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Abstract 

One  of  the  main  problem  affecting  modern  propulsor  design  engineers  is  the  ability  to 
quantitatively  predict  unsteady  propeller  forces  for  modern,  multi-blade  row,  ducted 
propulsors  operating  in  highly  contracting  fiowfields.  Current  algorithms  provide 
valuable  insight  into  qualitative  trendlines  for  these  modern  designs.  This  thesis 
has  focused  on  the  more  accurate  quantitative  force  prediction  by  introducing  more 
physical  modeling  into  the  numerical  computations,  using  more  accurate  analytical 
representation  of  continuous  physical  phenomena,  whilst  not  increasing  the  usage 
complexity  for  the  desktop  engineer. 

This  thesis  developed  several  novel  algorithms  and  techniques  and  applied  them 
to  build  an  evolutionary,  general  vortex-lattice  lifting-surface  propeller  code.  First, 
a  general  method  to  track  the  trajectory  of  individual  wake  singularity  sheets  and 
compute  their  influence  velocities  was  evolved  which  reduces  computation  time,  and 
dramatically  increases  the  accuracy  of  the  unsteady  blade  loading  problem. 

To  improve  the  general  coupling  technique  between  potential-based  propeller  codes 
and  volumetric  Reynolds-Averaged  Navier-Stokes  codes,  a  general  analytic  method 
based  upon  an  elliptic  integral  method  for  the  velocity  induced  by  a  vortex  ring  of 
unsteady  harmonic  strength  to  compute  of  the  time-averaged  induced  velocities  in  the 
swept  volume  of  the  propeller  was  introduced  which  is  more  accurate,  as  demonstrated 
in  model  problems,  and  more  robust,  as  indicated  by  improved  convergence  and  ac¬ 
curacy  in  a  fully  three  dimensional  propeller  code.  A  discretized  geometric  technique 
was  also  created  to  internalize  the  coupling  routines,  making  the  code  more  robust, 
while  decreasing  the  computation  burden  over  currect  methods. 

Finally,  a  higher  order  quadratic  influence  function  technique  was  implemented 
within  the  wake  to  more  accurately  define  the  induction  velocity  at  the  trailing  edge 
which  has  suffered  in  the  past  due  to  lack  of  discretization.  These  propeller  propgram 
enhancements  were  fitted  into  a  fully  functional  version  of  the  Propeller  Unsteady 
Forces  (PUF)-series  code,  and  coupled  with  a  three  dimensional  RANS  code. 

Thesis  Supervisor:  Justin  E.  Kerwin 
Title:  Professor  of  Naval  Architecture 
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Chapter  1 


Thesis  Overview  and  Motivation 


One  of  the  main  problem  affecting  modern  propulsor  design  engineers  is  the  ability  to 
quantitatively  predict  unsteady  propeller  forces  for  modern,  multi-blade  row,  ducted 
propulsors  operating  in  highly  contracting  flowfields.  Current  algorithms  provide 
valuable  insight  into  qualitative  trendlines  for  these  modern  designs.  This  thesis 
has  focused  on  the  more  accurate  quantitative  force  prediction  by  introducing  more 
physical  modeling  into  the  numerical  computations,  using  more  accurate  analytical 
representation  of  continuous  physical  phenomena,  whilst  not  increasing  the  usage 
complexity  for  the  desktop  engineer. 


1.1  Thesis  Goals 

The  goal  of  this  thesis  is  to  present  improved  algorithms  in  production  quality  software 
that  allows  naval  engineers  to  more  accurately  and  more  quickly  predict  the  unsteady 
forces  and  moments  on  advanced  propulsors.  Accomplishing  this  goal  will  allow  more 
flexibility  in  including  advanced  unsteady  analysis  as  routine  part  of  the  exploration 
of  the  design  space  presented  to  an  engineer  at  the  outset  of  any  project.  In  fact,  the 
bounds  upon  the  design  space,  at  the  outset,  are  more  determined  by  computation 
setup,  running  complexity  and  computing  time,  rather  than  other  inherent  challenges. 

The  goal  of  improved  algorithms  led  to  the  development  of  analytic  routines  to 
replace  discretized  routines  which  were  shown  to  have  slow  convergence  rates.  In 


17 


other  cases,  unstable  algorithms  were  replaced  with  more  stable  algorithms.  And 
physical  insights  allowed  for  major  speed-ups  with  no  loss  of  accuracy. 

The  goal  of  production  quality  software  means  that  internal  assumptions  must 
be  clearly  stated,  converged  and  validated,  or  offered  to  the  user  to  alter  in  an  easy- 
to-use  manner.  This  is  accomplished  in  this  thesis  through  the  use  of  input  files  for 
controlling  the  code,  vice  requiring  users  to  recompile  the  code  frequently. 

The  final  validation  of  these  goals  is  in  the  comparison  of  the  results  to  previous 
unsteady  propeller  codes,  and  experimental  results. 


1.2  Unsteady  Propeller  Forces 

The  accurate  prediction  of  the  steady  and  unsteady  propulsor  forces  arising  from 
operation  of  the  propulsor  in  non-uniform  wakes  is  an  important  aspect  in  the  pre¬ 
liminary  design  or  later  stage  analysis  of  any  propulsor.  By  unsteady  it  is  meant  that 
the  flow  appears  unsteady  to  a  bug  which  sits  upon  the  leading  edge  of  the  propeller. 
To  an  observer  in  the  global  reference  frame,  the  flow  will  appear  quite  steady,  but  the 
bug  will  see  large  changes  in  angles  of  attack  and  other  inflow  phenomena  traveling 
through  360  degrees. 

In  a  practical  sense,  the  predicted  unsteady  propeller  blade  forces  lead  to  numerous 
other  design  considerations: 

1.  Shaft  bearing  and  stern  tube  seal  requirements. 

2.  Determination  of  acceptable  vibration  and  radiated  acoustic  levels. 

3.  Steady  off-axis  maneuvering  forces  which  must  be  counteracted  by  other  control 
surfaces  and  drive  seal  and  stern  strength  considerations. 

4.  Unsteady  cavitation  inception,  which  again  affects  noise  and  blade  corrosion. 


18 


1.3  Historical  Context 


This  thesis  follows  from  and  complements  a  long  line  of  research  conducted  into  the 
prediction  of  propeller  forces  arising  during  operation  of  a  rotating  propeller  in  an 
unsteady  flowfield.  While  the  propeller  code  is  the  heart  of  this  thesis,  the  strong 
coupling  between  the  propeller  and  the  surrounding  flow  field  must  be  accounted  for 
by  the  use  of  an  external  flow  solver.  Thus,  the  issues  surrounding  accuracy,  ease 
of  use,  and  speed  of  use,  associated  with  coupling  the  propeller  and  flow  solvers  are 
addressed. 

The  propeller  code  used  as  a  starting  point  for  this  research  is  the  Propeller 
Unsteady  Forces  (PUF)  code  developed  at  the  Massachusetts  Institute  of  Technology 
(MIT)  originally  by  Kerwin  and  Lee  [16].  The  vortex-lattice  lifting-surface  theory 
which  forms  the  basis  for  this  research  is  based  up  the  theory  presented  by  Kerwin 
and  Lee  [16].  Recent  work  by  Warren  [30]  highlighted  the  need  for  improved  modeling 
of  the  wake  near  the  trailing  edge,  which  directly  led  to  the  work  presented  in  Section 
4.9.  Work  by  Kinnas  [17]  showed  that  for  contrived  inflows  wake  alignment  issues 
were  important  in  unsteady  blade  force  calculations. 

The  Reynold’s- Aver  aged  Navier-Stokes  code  used  for  this  research  is  IFLOW  de¬ 
veloped  at  David  Taylor  Model  Basin  by  Sung  [27], [28].  The  RANS  code  UNCLE, 
developed  at  Mississippi  State  University  was  also  used  to  verify  the  validity  of  the 
approaches  taken  here. 

1.4  Summary  of  Computational  Improvements 

This  thesis  has  resulted  in  significant  improvements  to  the  unsteady  analysis  of  propul- 
sors  by  providing  more  accurate,  more  computationally  efficient,  and  more  robust 
computer  analysis  programs. 
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transfer  PUF  forces  to  a  RANS  code 

Propeller  Unsteady  Forces  (PUF)  Program  and  RANS  Program  Coupling 


This  Research 

Previous  Methodology 

Analytic  method  for  prediction  of  induced 

velocity  via  ring  vortices  of  unsteady  har¬ 
monic  strength 

Discrete  algortihm  via  Biot-Savart  Law 

Fast  geometric  subdivision  algorithm  for 

transferring  forces 

Polygon  overlap  methodology 

Trilinear  interpolation  used  to  transfer  ve¬ 
locity  field  information  between  PUF  and 

RANS  domains 

Cubic  interpolation  scheme,  which  is  un¬ 
stable  in  regions  of  high  gradients 
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1.5  Summary  of  Results 


This  thesis  has  produced  an  improved  vortex-lattice  lifting-surface  analysis  code, 
PUF,  as  well  as  a  revised  RANS  code,  IFLOW,  which  ,  when  used  in  tandem,  give 
fast,  accurate  solutions  to  the  most  complex  of  geometries. 

•  The  new  PUF  is  a  vortex-lattice  lifting-surface  propeller  code  which  incorpo¬ 
rates  fully-three  dimensional,  independently  aligned  wakes. 

•  IFLOW  is  a  three  dimensional  RANS  code  which  internally  couples  with  a 
propeller  blade  force  solver. 

Both  of  these  codes  can  be  run  independently,  or  coupled  with  a  single  change  to 
an  input  file.  This  fiexibility  removes  several  degrees  of  freedom  from  the  analysis 
process  and  allows  engineers  to  concentrate  on  the  engineering  design  issues,  rather 
than  the  analysis  methodology  issues. 
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Chapter  2 


Vortex-Lattice  Lifting-Surface 
Formulation 


Vortex-lattice  lifting-surface  theory  models  a  propeller  as  a  singularity  distribution  in 
free  space,  and  imposes  the  kinematic  boundary  condition  at  selected  points  on  the 
foil  surface  requiring  there  to  be  no  flow  through  the  propeller  blade  at  that  point. 


2.1  Lifting  Line  Theory 

The  roots  of  vortex-lattice  lifting-surface  theory  lie  actually  in  the  development  of 
the  lifting  line  theory  early  in  the  20th  century  [5].  In  lifting  line  theory,  a  wing 
is  represented  as  a  two  dimensional  line  of  bound  vorticity,  as  shown  if  Figure  2.1. 
Figure  2.1  shows  a  wing  modeled  as  a  lifting  line  constant  strength  vortex.  Due  to 
Kelvin’s  thereom,  which  briefly  states  that  the  circulation  within  a  material  volume 
is  constant  over  time,  the  lifting  line  of  vorticity  which  lies  along  the  wing  must 
extend  downstream  of  the  wing  where  there  is  a  line  of  vorticity  of  equal  and  opposite 
magnitude  downstream  of  the  lifting  line  wing.  Taking  a  material  volume  as  any 
simple  enclosed  volume  which  surrounds  the  convex  hull  of  the  lifting  line  wing  and 
its  trailing  vortex  system  (composed  of  two  trailers  and  the  starting  vortex),  we  see 
that  Kelvin’s  thereom  is  satisfied  since  the  sum  of  all  the  circulation  due  to  the  lifting 
line  is  zero.  A  more  realistic  vortex-lattice  lifting  line  representation  of  the  wing  which 
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Figure  2-1:  Lifting  line  representation  of  a  finite  wing. 


has  non  constant  strength  along  the  span  can  be  implemented  by  stacking  individual 
vortex  segments  shown  in  Figure  2.1  end  on  end. 

A  mathematical  requirement  of  the  lifting  line  wing  representation  is  that  the 
circulation  on  the  tips  of  the  wings  be  zero.  A  completely  rigorous  mathematical 
derivation  is  presented  in  both  Newman  [22]  and  Kerwin  [14].  Due  to  this  singularity 
requirement,  the  circulation  distribution  over  the  span  of  a  lifting  line  is  usually 
represented  as  a  sine  series,  as  first  demonstrated  by  Glauert  [5]. 

The  lifting  line  wing  has  been  implemented  as  a  vortex-lattice  lifting-line  wing  by 
Kerwin  [14].  The  lifting  line  wing  theory  is  directly  transferable  to  propellers,  and 
was  implemented  as  the  Propeller  Lifting  Line  computer  program  by  Coney  [2], [3]. 
The  problem  with  the  lifting  line  representation  of  the  blade,  though,  is  that  it  can 
not  represent  swept  wings  (or  skewed  blades)  due  to  the  increasing  singularity  present 
in  the  Cauchy  principal  value  integral  equation.  The  lifting-surface  theory  removes 
this  singularity  by  spreading  the  concentrated  bound  vortex  over  a  chord. 

2.2  Vortex-Lattice  Theory 

In  a  traditional  propeller  lifting  surface  propeller  code,  a  grid  lattice  is  placed  on 
the  blade  mean  camber  surface,  the  hub  and  duct  (if  present)  and  the  trailing  wake 
system  [16].  Each  lattice  segment  is  assigned  a  strength  of  vorticity.  On  the  solid  body 
surfaces,  such  as  the  blade  and  hub,  a  control  point  is  placed  near  the  center  of  the 
grid  lattice.  The  strength  of  each  vortex  lattice  segment  is  assigned  by  satisfying  the 
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kinematic  boundary  condition  that  the  flow  must  be  tangent  to  the  propeller,  hub,  and 
duct  surfaces  at  every  control  point.  A  wake  lattice  is  traditionally  found  by  laying 
the  lattice  upon  a  helicoidical  surface,  where  the  pitch  is  set  by  the  hydrodynamic 
pitch  possibly  with  induced  effect  corrections. 

Mathematically,  the  propeller  problem  involves  a  simple  matrix  equation.  By 
attacking  the  geometry  of  the  problem,  an  influence  matrix  is  formulated  which  gives 
the  velocity  induced  by  a  unit  strength  of  vorticity  along  every  vortex  lattice  upon 
every  control  point  -  an  [/ATFjluence  matrix. 


[INF] 


—  ^induced 


(2.1) 


The  symmetry  of  the  blades  and  wakes  allows  for  great  computational  simplification 
through  similarity  principles.  Next,  to  satisfy  the  kinematic  boundary  condition,  the 
physics  of  the  problem  dictate  that  the  component  of  the  induced  velocity  normal  to 
the  blade,  when  added  to  the  component  of  the  effective  inflow  velocity  normal  to  the 
blade  must  be  zero. 


+  Ve 


eff  •  — 


(2.2) 


Equation  (2.2)  is  the  heart  of  the  propeller  lifting  surface  code. 

Equation  (2.2)  highlights  the  fact  that  a  valid  propeller  calculation  requires  both 
the  correct  geometrical  position  of  the  lattice  and  the  correct  effective  velocity  every¬ 
where  on  the  blade  surface. 

The  propeller  forces  resulting  from  the  vorticity  and  source  distributions  are  calcu¬ 
lated  from  the  Kutta-Joukowski  and  Lagally  theorems,  respectively.  Unsteady  forces 
arise  from  the  unsteady  pressure  term  in  the  unsteady  Bernoulli  equation.  A  Lighthill 
leading  edge  suction  force  correction  is  applied  to  these  forces  in  a  chordwise  fash¬ 
ion,  and  the  propeller’s  sectional  drag  is  calculated  either  based  on  strip  wise  two 
dimensional  empirical  drag  coefficients,  or  a  stripwise  2-D  integral  boundary  layer 
calculation  [9].  Cavitation  effects  can  be  included  through  the  use  of  a  pressure  and 
velocity  indication  which  finds  the  detachment  and  re-attachment  points  of  the  cavity 
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Lattice  of 
Bound  Vortex 


Figure  2-2:  A  vortex-lattice  lifting-surface  propeller  blade  and  wake  model. 


sheet.  The  computer  code  developed  in  this  thesis  uses  only  empirical  sectional  drag 
coefficients  and  does  not  account  for  cavitation  effects. 
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Chapter  3 


Unsteady  Vortex-Lattice 
Lifting-Surface  Formulation 


In  an  unsteady,  vortex-lattice  lifting-surface  propeller  code  the  requirement  to  satisfy 
both  a  Kutta  condition  and  Kelvin’s  thereom  drives  scientists  to  devise  ever  new 
variations  on  a  numerical  theme  to  align  theory  with  observation. 

The  main  feature  of  the  lifting  surface  formulation  requires  a  representation  of  a 
lifting  body. 

•  Camber  effects  are  accounted  for  by  modeling  the  lifting  surface  on  the  mean 
camber  surface. 

•  Induced  effects  due  to  both  bound  and  shed  vorticity  are  calculated. 

•  Angle  of  attack  effects  are  found  by  properly  modeling  the  three  dimensional 
surface  rotating  through  the  fluid. 

The  paper  thin  vortex  lattice  representing  the  blade  is  set  upon  the  mean  camber 
surface.  Consequently  provisions  must  be  made  for  the  two-dimensional  effects  of  the 
foil: 

•  Thickness.  Can  be  accounted  for  through  the  distribution  of  source  singularities 
alongside  the  vortex  singularities. 
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•  Leading  edge  suction.  Foils  not  at  their  shock  free  entry  angle  produce  leading 
edge  suction.  This  can  be  handled  through  the  use  of  an  analytically  correct 
2-D  solution,  applied  in  three  dimensions,  as  shown  by  Lan  [19]. 

•  Boundary  layer  (shear)  drag.  Any  number,  or  combination  of,  empirical,  semi- 
empirical,  or  2-D/3-D  analytic  methods  can  be  implemented  [9]. 

3.1  Kutta  Condition 

In  general,  the  Kutta  condition  can  take  many  forms: 

1.  Zero  pressure  jump  at  the  trailing  edge. 

2.  Finite  velocity  at  the  trailing  edge. 

3.  Smooth  flow  tangent  to  the  blade  at  the  trailing  edge  [22]. 

4.  A  wake  singularity  sheet  extends  downstream  of  the  blade  where  the  singularity 
strength  is  the  difference  in  singularity  strengths  on  the  upper  and  lower  surfaces 
of  the  blade  directly  upstream  of  the  trailing  edge 

The  numerical  implementation  of  the  Kutta  condition  then  falls  classically  into  either 
an  explicit  or  an  implicit  form  [16]. 

In  the  explicit  formulation  of  the  Kutta  condition  the  strength  of  the  last  few 
bound  vortices  on  the  blade  are  manipulated  to  satisfy  a  condition  imposed  by  the 
researcher.  For  steady  flow,  the  condition  that  the  bound  vortex  strength  vanish 
at  the  trailing  edge  is  imposed.  Numerically,  the  imposition  of  this  extra  condition 
requires  either  another  unknown,  or  that  one  of  the  constraint  equations  be  dropped. 

In  the  implicit  formulation  of  the  Kutta  condition,  judicious  placement  of  the 
blade  and  wake  vortices  are  assumed  to  satisfy  the  Kutta  condition.  The  hypothesis 
can  be  assured  by  placing  the  final  downstream  control  point  along  the  trailing  edge. 
Since  the  numerics  of  the  discretized  boundary  value  problem  require  there  to  be  no 
fluid  velocity  normal  to  the  blade  at  any  control  point,  the  Kutta  condition  is  assured. 
Further  derivations  and  supporting  arguments  have  been  shown  in  the  literature. 
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Steady  Formulation 

Free  Vorticity  -  Horseshoe  extends  to  infinity 


7o  Jo 

F  Starting  vortex 

^  is  infmtely  downstream 

Unsteady  Formulation 

Free  Vorticity  -  Loops  correspond  to  time  step 

/9°  /° 

_2!_ 

/ir.-  /ir- 

//  *  / 

//  //  / 

//  y  / 

//  //  / 

Vortex  Loop  Spacing  on  Blade  set  by  user  Vortex  Loop  Size  in  Wake  depends  on  Total  Flow  Velocities 

Figure  3-1:  These  two  figures  show  the  major  difference  between  steady  (upper  figure) 
and  unsteady  (lower  figure)  vortex  lattices  formulations. 

3.2  Unsteady  Modeling 

The  major  difference  between  a  steady  and  unsteady  formulation  is  the  shedding  of 
span-wise  vortex  segments  of  free  vorticity  which  account  for  the  changes  in  blade 
circulation  as  the  propeller  traverses  three  hundred  and  sixty  degrees.  This  is  shown 
in  Figure  3-1.  In  a  steady  code  the  wake  vortex  is  a  single  horseshoe  extending  to 
infinity.  In  an  unsteady  code  the  wake  vortex  loop  strength  changes  due  to  the  change 
in  blade  shed  vorticity  at  each  blade  time  step.  So  along  a  constant  spanwise  series  of 
wake  vortex  loops,  there  are  unequal  amounts  of  vorticity.  The  difference  in  strength 
of  the  trailing  vortex  segments  (*.e.,  the  vortex  lattice  segments  in  the  streamwise 
direction)  require  that  spanwise  vortex  segments,  called  shed  vortex  segments,  be 
introduced  which  connect  the  trailers  to  satisfy  Kelvin. 

In  the  unsteady  code,  the  unknowns  include  not  only  the  vortex  strength  on  the 
blades,  but  also  the  first  wake  loop.  This  is  because  as  the  blade  advances  through 
the  fluid  from  time  step  N  to  time  step  iV  -I- 1  the  amount  of  circulation  shed  into 
the  wake  to  satisfy  Kelvin  is  unknown.  So  the  general  solution  procedure,  as  shown 
in  Figure  3-2  is 

1.  Starting  from  a  known  solution,  advance  the  wakes  downstream  by  one  loop. 
Note  that  we  make  the  assumption  that  the  loops  are  sized  such  that  one  loop  is 
equivalent  in  size  to  the  distance  a  fluid  particle  would  travel  in  one  time  step. 
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Figure  3-2:  This  figure  shows  a  slice  along  the  chord  wise  direction  (constant  span) 
for  a  blade,  and  the  steps  required  to  advance  the  blade  one  time  step. 


2.  Solve  the  boundary  value  problem  on  the  blade  surface. 

3.  Place  the  correct  circulation  in  the  first  wake  loop  based  on  the  definition  of 
Kelvin’s  thereom. 


4.  This  then  represents  another  known  solution,  and  the  process  repeats  itself. 

One  requirement  to  implement  this  scheme,  is  that  the  vortex  lattice  spacing  in  the 
wake  is  directly  dependent  upon  the  time  step  set  by  the  user.  The  blade  vortex 
spacing,  though,  is  set  in  a  different  manner  by  the  user  since  there  is  no  intra¬ 
problem  dependency  between  the  number  and  relative  size  of  the  blade  and  wake 
vortex  lattices.  This  gives  rise  to  a  sharp  gradient  in  grid  density  at  the  trailing 
edge  of  the  blade  as  the  solution  transitions  from  the  blade  to  the  wake  lattice.  This 
problem  is  further  accentuated  by  the  fact  that  the  blade  lattice  is  cosine  spaced  in 
the  chordwise  direction. 

staS^tlitm 

vortex  induced  flow  due  the  vorticity  on  the  blade,  the  vortex  induced 
flow  due  to  the  vorticity  in  the  wake,  and  the  effective  inflow. 

3.3  Discrete  Bjpundary  ,yalue  Problem 

=  -V,  (3.1) 

j=l  k=l 
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where 


Nb  number  of  blade  vortex  loops 
(Fft)  j  circulation  strength  of  blade  vortex  loop  j 
Bij  Infuence  of  unit  circulation  strength  in  blade  j  on 
control  point  i 

Ny;  number  of  wake  vortex  loops 
(rs)fc  circulatoin  strength  of  wake  loop  k 
Wi^k  Infuence  of  unit  circulation  strength  in  wake  loop 
k  on  control  point  i 

Vi  Effective  velocity  at  control  point  i 


In  the  unsteady  formulation,  a  simple  notation  scheme  is  introduced  modifying  Equa¬ 
tion  3.1  which  states  that  the  end  state  of  the  problem  is  to  find  the  solution  at  the 
next  time  step,  -t- 1. 


Nb  iVw 

53(r,)''+‘By  +  =  -K,  (3.2) 

j=l  k=l 

From  the  physics  of  the  problem,  it  is  known  that  the  vorticity  present  in  the  wake 
convects  with  the  local  velocity.  By  sizing  the  wake  panels  such  that  the  size  of  one 
wake  panel  is  the  distance  traveled  by  the  blade  in  one  time  step,  we  can  state  with 
certainty  that  the  vorticity  present  on  the  A:  +  1  panel  at  time  iV  -f  1  is  the  vorticity 
which  existed  on  the  k*"^  panel  at  time  N. 


(r.)K' =  (r.)J'  for  A:>0  (3.3) 

Since  the  shed  vorticity  present  in  the  wake  is  known  at  time  N,  this  means  that  the 
wake  influence  velocities  are  known  due  to  all  wake  panels  except  the  k  =  1  panel. 


Nt 


Nu, 


+  (r,)^'^^  =  -K  -  5;(r,)f 


(3.4) 


j=l 


fc=2 
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To  account  for  the  vorticity  present  on  the  first  wake  panel, we  make  use  of 
Kelvin’s  thereom.  When  applied  to  a  material  contour  which  encompasses  the  blade 
and  the  entire  wake,  we  can  see  that  any  change  in  vorticity  on  the  blade  over  a  given 
time  interval  must  be  balanced  by  a  change  in  vorticity  in  the  wake. 


+ 


dt 

5r, 


dt  dt 


0 

0 


(3.5) 


The  first  observation  to  make  concerning  Equation  3.5  is  that  integrating  all  terms 
with  respect  to  time  t  gives  the  absolute  change  in  vorticity.  That  being  said,  the 
second  observation  to  make  is  that  the  only  place  the  vorticity  in  the  wake  changes 
within  the  material  volume  (which  by  the  definition  of  a  material  volume  convects 
with  the  flow)  is  on  the  first  wake  panel. 


dFf, 

dt 


dt  +  (Fj) 


N+l 

1 


=  0 


(3.6) 


Applying  a  very  low  order  integration  scheme,  such  as  a  one-sided  trapezoidal  rule 
for  the  first  term  in  Equation  3.6 


L  + 

Tft  =Er6 

ATfe 

pN+1  op 

=  ECr.)"*'  (3.7) 

Substituting  Equation  3.7  into  Equation  3.6  gives 
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•p  A^+1  _ 

-L  SI  — 


(3.8) 


iv. 

-E(r»)r' 

J=1 

Substituting  Equation  3.8  into  Equation  3.4 


Nb  Nb  Ny, 

(3-9) 

j=l  j=l  k=2 

Equation  3.9  is  the  final  discretized  boundary  value  problem  which  can  be  solved 
numerically.  The  important  things  to  note  about  Equation  3.9  is  that  all  the  un¬ 
knowns  on  the  left  hand  side  are  the  Ft  of  interest  at  the  next  time  step.  Further, 
the  low  order  integration  means  that  interestingly  enough  that  first  wake  panel  does 
not  influence  the  right  hand  side  of  the  equation. 

Combining  the  two  terms  on  the  left  hand  side  of  Equation  3.9  reduces  it  to  the 
more  usual  form  of  Ax  =  b. 


Nb 

E(r6)f^MBi3  -  M'i.il  =  -v)  -  (a-ic) 

j=l  k=2 

Physically  this  means  that  the  blade  influence  loops  include  the  influence  of  the  first 
wake  panel.  Note  that  the  influence  of  the  first  wake  panel  is  not  found  on  the  right 
hand  side  of  Equation  3.10. 

This  formulation  differs  from  previous  unsteady  formulations  in  the  placement  of 
the  first  shed  vortex  segment.  Previous  researchers  have  placed  the  first  shed  vortex 
segment  one  quarter  downstream  of  the  foil  trailing  edge,  presumably  in  an  attempt  to 
avoid  the  trailing  edge  singularity.  This  thesis  has  shown  that  the  correct  placement 
of  the  first  shed  vortex  is  one  complete  time  step  downstream  of  the  blade  trailing 
edge. 
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3.4  Higher  Order  Wake  Representation 


One  method  to  increase  the  accuracy  for  a  given  problem  is  to  increase  the  order  of  the 
discretized  representation  of  the  physical,  analog,  system.  For  a  vortex-lattice  code, 
this  strategy  leads  to  a  higher  order  representation  of  the  singularity  distribution  over 
a  given  vortex  loop  since  in  its  most  basic  form,  a  vortex  loop  is  a  constant  strength 
dipole  sheet.  Practically  speaking,  this  is  implemented  by  subdividing  a  vortex  loop 
and  assigning  a  spatially  varying  singularity  strength  to  the  subdivided  loops  based 
on  the  singularity  value  at  the  boundary  of  the  original  loop. 

For  a  lifting  surface  code,  there  is  no  requirement  to  represent  the  entire  transition 
wake  as  higher  order  loops,  because  most  of  the  loops  are  far  downstream  from  the 
control  points  of  interest.  To  save  computational  time,  then,  only  the  first  wake 
panel  downstream  of  the  trailing  edge  is  represented  by  a  higher  order  singularity 
distribution.  This  solution  technique  assumes  that  the  blade  solution  is  affected  most 
by  the  order  of  the  first  wake  panel  since  the  first  wake  panels  adjacent  to  the  trailing 
edge  have  the  largest  influence  upon  induced  velocities  at  the  boundary  value  problem 
control  points. 

The  theory  and  results  are  presented  below  for  a  linear  variation,  and  2nd  order 
variation  in  singularity  strength. 


3.4.1  Linear  Variation  in  Wake  Loop  Strength 

The  first  improvement  upon  the  low  order  scheme  presented  in  Equation  3.10  is  to 
model  the  strength  of  vorticity  across  the  first  wake  panel  as  linearly  varying  across 
the  panel. This  is  accomplished  by  breaking  the  first  shed  wake  panel  into  several 
smaller  panels  as  shown  in  Figure  3-3.  Higher  order  wake  panels  are  modeled  by 
subdividing  the  first  wake  panel.  The  singularity  strength  of  each  sub  panel  can  then 
be  forced  to  conform  to  a  higher  order  scheme:  linear,  quadratic  or  higher.  Notice 
that  the  strength  across  the  panel  then  depends  upon  the  blade  circulation  which  is 
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First  wake  loop  subdivided  into  m  subpanels 

Figure  3-3:  Higher  order  wake  panel. 

the  unknown,  and  the  previous  timestep  blade  circulation,  which  is  known. 

=  (1 -m)r,o  +  mr„ 

r.r'  =  (1  -  m)  Efi  +■  +  (m)  Efi  (xii) 

The  variable  m  in  Equation  3.11  represents  the  distance  between  0  and  1  across  the 
subdivided  panel. 

This  changes  the  discretized  boundary  value  problem  as  shown  in  Equation  3.12. 
Notice  in  Equation  3.12  that  the  influence  of  the  first  wake  panel  is  accounted  for  on 
both  the  left  hand  side  and  right  hand  side  of  the  equation.  And  that  each  subdivided 
panel  influences  each  control  point. 

Ni,  Nm  Nu,  Nm  JV(, 

j-l  m-1  k=2  m=l  j=l 

(3.12) 

where  Wm,i  is  the  influence  of  the  subpanel  in  the  wake  panel 

3.4.2  Quadratic  Variation  in  Wake  Loop  Strength 

The  second  improvement  upon  the  low  order  scheme  presented  in  Equation  3.10  is  to 
model  the  strength  of  vorticity  across  the  first  wake  panel  as  quadratically  varying. 
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This  is  accomplished  by  again  breaking  the  first  shed  wake  panel  into  several  smaller 
panels  as  shown  in  Figure  3-3. 

Examining  a  single  chordwise  strip  of  vortex  loops  extending  from  the  blade  and 
into  the  wake,  this  method  states  that  the  strength  of  vorticity  in  the  wake  panel 
varies  quadratically  in  space  downstream  of  the  blade  trailing  edge.  Therefore,  the 
strength  of  the  vorticity  is  represented  as  a  quadratic  equation  as  shown  in  Equation 
3.13. 


r "h  bs  -|-  c  (3.13) 

where  s  =  distance  downstream  of  trailing  edge 
a,b,c  =  polynomial  coefficients 

Because  the  wake  panel  spacing  is  constant,  a  3X3  matrix  can  be  formed  to  find 
the  polynomial  coefficients  using  the  unknown  T^o  at  the  trailing  edge,  and  the  known 
Tsi  and  rs2  on  the  second  and  third  panels.  When  the  coefficient  matrix  is  inverted 
and  the  terms  re-ordered  to  show  the  dependence  on  the  various  T^’s,  Equation  3.14 
results.  This  is  equivalent  to  solving  the  Lagrange  polynomial  solution. 


By  following  up  with  the  usual  representation  of  the  shed  vorticity  in  terms  of  the 
blade  vorticity  at  a  given  time  step 


Ni 

{^s)te  =  y^(rt) 

Nb 

r.i  =  £(r.)' 

j~l 


iV+1 

3 


N 


36 


Nb 

V.. = E(r.)r- 

j=l 

Inserting  these  definitions  into  Equation  3.10  and  migrating  the  unknowns  to  the 
left  hand  side  and  the  known  to  the  right  hand  side  leads  to  the  final  discrete  equation 
shown  in  Equation  3.15. 


Nb  Nm  ^  „ 

--S  + = 


j=l 


m=l 
Nm  Nb 


Nb 


k=2 


m=l  j=l 


(3.15) 

where  s  is  the  non-dimensional  distance  of  the  sub¬ 
divided  wake  panel  across  the  total  wake 
panel 
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Chapter  4 


Force  Free  Wake  Alignment 


Force  free  wake  alignment  means  that  the  wake  vortex  lattice  structures  are  aligned 
with  the  local  fluid  velocity.  The  wake  is  then  force  free  by  the  Kutta-Jukowski 
thereom,  Equation  4.1,  since  the  cross  product  of  the  velocity  and  vorticity  is  zero. 

'f  =  pVx^  (4.1) 

This  chapter  presents  the  development,  implementation,  and  validation  of  a  novel 
technique  to  align  the  vortex-lattice  wake  with  the  local  flowfield.  The  engineering 
applicability  to  this  line  of  research  is  presented  in  Section  4.1.  A  broad  overview  to 
wake  modeling  in  the  context  of  a  propeller  code  is  presented  in  Section  4.2.  Section 
4.3  presents  the  approach  taken  to  accomplish  three  dimensional,  independent  wake 
alignment.  Two  model  problems  are  shown  in  Section  4.4.  Section  4.5  contains  the 
details  on  wake  growing.  Hub  image  placement  is  shown  in  Section  4.6.  A  very 
relevant  adjunct  to  wake  alignment  which  makes  the  computational  burden  bearable 
is  wake  subdivision,  shown  in  Section  4.8.  A  novel  method  to  improve  computation 
accuracy  at  reduced  computational  expense  through  judicious  application  of  higher 
order  techniques  is  presented  in  Section  4.9,  where  a  higher  order  wake  panel  to  resolve 
the  influence  velocity  more  accurately  near  the  trailing  edge  singularity  is  created. 
Results  for  these  new  algorithms  as  implemented  in  a  fully  functional  propeller  code 
are  found  in  Section  4.10. 


39 


4.1  Applications  for  Non-Uniform  Wakes 


Real  world  applications  of  propeller  codes  highlight  two  cases  of  interest  where  wake 
alignment  contributes  to  the  accuracy  of  the  solution: 

1.  Propellers  operating  on  inclined  shafts  and  the  non-axial  flow  which  follows  the 
aft  buttock  lines. 

2.  Maneuvering  vehicles  present  strong  cross  flow  contributions  to  the  inflow  seen 
by  the  propeller  in  its  rotating  reference  frame. 

3.  Propellers  operating  within  unsteady  flowfields  due  to  the  presence  of  upstream 
appendages  such  as  control  surfaces,  sails,  struts,  or  an  inclined  shaft. 

These  real  world  flow  conditions  lead  to  unsteady  propeller  blade  forces.  The 
naval  engineer  then  has  the  task  of  designing  the  propeller  such  that  the  outwardly 
measurable  effects  of  the  generated  unsteady  blade  forces  can  be  accurately  modeled, 
in  magnitude  and  variation,  to  within  the  design  specifications. 

There  are  two  approaches  that  can  be  taken  to  align  the  propeller  vortex-lattice 
wake  with  the  local  fluid  velocities: 

1.  Align  the  vortex  lattice  wake  using  the  combined  affects  of  the  blade  and  wake 
singularity-induced  velocity  and  a  given  effective  inflow. 

2.  Align  the  vortex  lattice  wake  with  a  prescribed  total  inflow.  The  total  inflow 
can  be  either  recorded  in  an  experimental  test  facility  or  found  with  a  numerical 
simulation  technique  such  as  the  RANS  used  in  this  research. 

The  first  approach  has  been  undertaken  by  Shih  [25]  and  Kinnas  [23]  among  others. 
This  approach  brings  with  it  the  requirements  to  lump  vorticity  at  the  tip  of  the 
wake  sheet  as  the  tips  roll-up  and  to  de-singularize  the  influence  kernel  function.  This 
approach  also  has  the  more  serious  shortcoming  that  it  does  not  properly  account  for 
the  fact  that  the  vortex-vortex  interaction  between  the  vorticity  present  in  the  inflow 
and  the  blade  vorticity  alters  the  inflow  velocity  profile  [24].  This  is  explained  more 
fully  in  Appendix  B.l. 
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The  second  approach,  that  of  aligning  the  the  vortex  lattice  wake  with  a  prescribed 
velocity  field  is  undertaken  in  this  thesis.  This  research  utilizes  the  basic  model  of 
coupling  with  a  Reynolds- Averaged  Navier-Stokes  flow  solver  [15]  to  give  directly  the 
wake  field.  There  are  two  main  reasons  as  to  the  relevance  of  this  approach: 

1.  The  general  fidelity,  usability,  and  availability  of  RANS  codes  makes  them  ac¬ 
cessible  to  design  engineers. 

2.  A  RANS  code  can  properly  model  the  vortex-vortex  interactivity,  redistributing 
the  vorticity  present  in  the  inflow  in  the  presence  of  the  propeller. 

One  possible  shortcoming,  which  is  not  explored  further,  is  that  the  time-averaged 
nature  of  the  RANS  solution  means  that  the  fluid  velocities  are  time- averaged  (  or 
equivalently  spatially-averaged),  while  the  lifting-surface  formulation  requires  local 
velocities  to  correctly  solve  the  boundary  value  problem  [15].  The  details  are  explained 
in  Appendix  B.  This  assumption  has  proven  valid  in  practice  for  axisymmetric  flow, 
where  circumferential  mean  flow  averaged  over  the  entire  propeller  swept  volume,  is 
substituted  for  local  blade  flow. 

There  are  two  major  advantages  of  the  approach  taken  here: 

1.  Computational  efficiency:  The  calculation  of  the  blade  and  wake  singularity 
influence  on  the  wake  sheets  is  an  extremely  time  consuming  process  which 
scales  with  the  square  of  the  number  of  singularity  elements.  Memory  can 
not  be  traded  against  the  computation,  since  the  influences  change  at  each 
alignment  step. 

2.  Compatibility  with  empirical  measurements:  Empirical  measurements  of 
flow  around  bodies  in  the  plane  of  the  propeller  provide  either  nominal  (velocity 
measurements  in  the  absence  of  the  propeller)  or  total  (velocity  measurements 
with  the  effects  of  the  propeller)  velocities.  For  use  in  numerical  validation,  if 
total  velocities  are  measured  these  can  be  exactly  used  as  boundary  conditions. 
Nominal  values  measured  in  an  empirical  test  can  be  used  to  calibrate,  or  set  the 
boundary  conditions,  for  a  numerical  solver.  In  both  cases,  if  the  measurements 
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are  taken  far  enough  upstream  of  the  propeller  to  negate  the  influence  of  the 
propeller  (usually  2  propeller  diameters  for  open  propellers)  these  measurements 
are  equal  and  can  be  used  as  a  boundary  condition  on  a  volume  solver.  Note, 
however,  that  this  method  assumes  that  the  propeller  code  is  coupled  with  an 
external  code  which  solves  for  the  effective  inflow  while  maintaining  the  proper 
velocities  downstream  of  the  propeller  to  convect  the  wake. 


4.2  Wake  Modeling  in  a  Propeller  Code 

The  wake  scheme  used  as  a  starting  point,  is  that  devised  by  Greeley  and  Kerwin  [7] 
as  implemented  in  both  the  PBD-14  code  [21]  and  the  original  PUF  code  [13].  In  this 
formulation,  the  wake  is  divided  into  both  a  transition  and  ultimate  wake  model. 

The  transition  wake  is  a  fully  three  dimensional  wake  vortex  lattice.  The  lattice  of 
the  transition  wake  matches  the  spanwise  number  of  blade  lattice  element,  while  the 
streamwise  number  of  vortex  loops  depends  upon  the  time  step  set  for  the  particular 
problem.  The  transition  wake  extends  downstream  from  the  blade  a  user-specified 
distance. 

The  ultimate  wake  model  is  composed  of  trailing  helical  vortices  extending  to 
infinity.  That  is,  at  each  spanwise  lattice  position,  a  helical  vortex  is  grown,  with 
the  same  pitch  as  the  tipwise  element.  The  closed  form  solutions  due  to  Hough 
and  Ordway  [8]  as  implemented  by  Leibman  [20]  provide  accurate  influence  velocities 
due  to  the  ultimate  wake  helical  vortices.  The  helical  vortices  at  a  given  radius  are 
averaged  at  each  blade  position,  to  produce  a  steady  set  of  helical  vortices. 

4.3  Approach  to  Wake  Alignment 

Two  pieces  are  required  to  solve  the  wake  alignment  problem: 

•  individual  wake  growing. 

•  faster  numerical  algorithms  to  compensate  for  the  increased  number  of  influence 
calculations  required. 
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The  first  piece  is  the  need  to  align  each  wake  vortex-lattice  sheet  convecting  down¬ 
stream  from  each  blade  with  the  local  flowfield.  Because  the  symmetry  of  the  axisym- 
metric  wake  generation  scheme  is  destroyed,  there  is  an  increase  in  computational  and 
memory  requirements  (since  each  wake  must  be  individually  grown  and  has  individ¬ 
ual  influence  calculations  and  associated  storage  requirements).  This  leads  to  the 
requirement  to  divorce  wake  discretization  from  problem  size.  This  is  accomplished 
through  the  implementation  of  a  pseudo-higher  level  wake  discretization  scheme  using 
sub-divided  wake  vortex  loops. 


4.4  Wake  Alignment  Model  Problems 


The  first  problem  for  the  two  cases  presented  in  Section  4.1  involves  the  propeller 
acting  in  the  non-uniform  flowfield  produced  by  upstream  appendages.  As  is  clearly 
shown  in  the  extreme  case  of  a  step  change  in  inflow  velocity  in  Figure  4-1  the  pitch 
angle  of  the  wake  is  directly  aflFected  by  the  non-uniformity  in  axial  velocity.  To 
quantify  these  effects  upon  the  resulting  propeller  higher  harmonic  forces,  a  prescribed 
inflow  composed  of  known  harmonic  inflows  can  be  imposed,  and  then  the  resulting 
harmonic  forces  predicted  by  the  propeller  program  with  aligned  and  non-aligned 
wakes  can  be  measured  and  compared. 

The  second  problem  of  wakes  convecting  downstream  at  angle,  as  in  the  case  of 
a  maneuvering  vehicle  or  a  propeller  mounted  at  a  fixed  angle  of  inclination  with 
respect  to  the  inflow  as  shown  in  Figure  4-2  ,  again  results  in  a  change  to  the  wake 
pitch  angle.  The  once  per  revolution  change  in  pitch  directly  affects  the  blade  solution. 
This  change  results  from  the  azimuthal  component  of  the  velocity  field  convecting  the 
wake.  Conceptually,  if  the  cartesian  representation  of  any  given  velocity  is  translated 
to  a  cylindrical  form,  and  radial  velocities  would  of  course  not  affect  the  pitch,  but 
axial  and  azimuthal  velocity  components  would  change  the  pitch. 
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Top  Half  Plane 
Velocity  =  1.5 


High  Pitch 


Lower  Half  Plane 
Velocity  =  0.5 


Low  Pitch 


Figure  4-1:  A  non-uniform  inflow  directly  affects  the  pitch  angle  of  the  wake  sheets. 
Here  the  non-uniformity  is  produced  by  a  step  function  change  in  the  axial  velocity. 


Figure  4-2:  With  the  propeller  operating  in  a  10  degree  inclined  inflow,  the  wake 
properly  follows  the  inflow  angle. 
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Figure  4-3:  As  the  wake  convects  downstream  aligned  with  the  local  flowfield,  it  forms 
a  cylindrical  tube. 

4.5  Wake  Tracking  Methodology 

The  PUF  reference  frame  is  stationary,  with  the  axis  centered  near  the  blade,  at 
the  centroid  of  blade  rotation.  Because  the  blades  remain  stationary  a  rotational 
component  is  added  to  the  transition  wake  as  it  convects  downstream  from  the  blade 
trailing  edge.  The  approach  taken  here  is  to  convect  the  wake  sheets  with  the  local 
flowfield.  As  the  wake  sheets  convect,  a  center  of  rotation  is  defined  at  each  timestep. 
The  wake  sheets  can  then  be  rotated  about  that  axis  by  the  angular  blade  step. 

If  the  wake  sheets  can  be  considered  to  convect  downstream  within  a  cylindrical 
tube  that  is  aligned  with  the  wakefield,  then  there  must  be  a  centerline  axis  along 
which  the  cylinder  can  be  generated.  This  is  shown  in  Figure  4-3.  Conversely,  if  the 
centerline  axis  (generator  line)  of  the  cylindrical  tube  can  be  defined,  then  the  wake 
sheets  can  be  properly  rotated  about  the  generator  line  to  account  for  the  rotational 
component  required  within  the  PUF  reference  frame.  Notice  for  a  propeller  in  a 
uniform,  straight  ahead  flow  situation,  like  in  a  water  tunnel,  the  wake  generator  line 
lies  along  the  shaft  axis.  In  real  world  cases  where  the  wake  shifts  to  align  with  an 
unsteady,  non-uniform  flow  field,  the  wake  generator  line  is  arbitrarily  oriented  in 
space.  This  leads  to  the  requirement  to  rotate  a  point  in  space  a  given  angle  about 
an  arbitrarily-oriented  axis. 
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Figure  4-4:  Notation  used  in  specifying  how  a  point  in  the  PUF  reference  frame  is 
rotated  about  an  arbitrary  vector 

4.5.1  Wake  Sheet  Generator  Line 

The  wake  sheet  generator  line  is  created  by  tracking  a  cylindrical  shell  of  fluid  through 
space.  The  centerline  of  the  cylindrical  shell  is  the  wake  sheet  generator  line.  The 
cylindrical  shell  is  generated  by  marching  a  ring  downstream  from  the  seven-tenths 
radius  of  the  propeller.  The  center  of  this  ring  is  the  centerline  axis.  The  position  of 
the  generator  line  at  each  timestep  is  found  by  averaging  the  coordinates  of  the  ring. 
The  orientation  of  the  generator  line  is  found  using  a  backwards  Euler  differencing 
scheme. 

Once  the  wake  sheet  generator  line  is  deflned,  the  wake  sheets  are  convected  with 
the  local  flow  and  rotated  by  the  propeller  angular  increment,  which  is  dependent  on 
the  number  of  propeller  time  steps  per  revolution.  This  involves  a  three  dimensional 
rotation  of  a  point  about  an  arbitrary  axis. 


4.5.2  Rotation  of  a  Point  about  an  Arbitrary  Axis 

The  rotation  of  a  point  {x,y,z)  about  an  arbitrary  vector  {ai,bj,ck}  positioned  in 
space  at  the  point  (xp,  yp,  zp)  by  a  desired  angle  a  is  found  using  a  series  of  rotational 
and  translational  matrix  manipulations.  The  basic  series  of  steps,  and  associated 
transformation  matrices,  are  shown  below: 
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•  Translate  the  coordinate  center  of  the  axis  of  rotation  to  the  axes  center. 


1  0  0  0 
0  10  0 

T  = 

0  0  10 

-XP  -YP  -ZP  1 

•  Rotate  the  axis  of  rotation  vector  about  the  Y-axis  to  place  the  coordinate 
center  of  the  vector  in  the  YZ  plane 

cos{9)  0  sin{9)  0 
0  10  0 

Re  = 

—sin{9)  0  cos{9)  0 
0  0  0  1 

•  Rotate  the  axis  vector  about  the  X-axis  to  align  the  rotation  axis  vector  with 
the  Z-axis. 

10  0  0 
0  cos{(^)  sin{4>)  0 

P'<j)  ~~ 

0  —sm{(f))  cos{4>)  0 
0  0  0  1 

•  The  rotation  axes  vector  and  the  global  Z  axes  are  now  colinear.  The  desired 
rotation  is  now  a  rotation  of  the  point  about  the  Z-axis  by  the  desired  angle  a. 

cos(ol)  sin{a)  0  0 
—sin(a)  cos(a)  0  0 
0  0  10 

0  0  0  1 

The  remaining  steps  remove  the  transformations  accomplished  in  steps  1-3. 
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•  Undo  the  rotation  about  the  X-axis 


R-6  — 


•  Undo  the  rotation  about  the  Y-axis 


R-fi 
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0 

0 

1 
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•  Translate  the  rotation  vector  base  point  (which  really  defines  the  local  coordi¬ 
nate  axes)  from  the  global  axes  center  to  original  position 

0  0  0 

1  0  0 

0  1  0 

[XP  YP  ZP  1 

These  matrices  can  be  combined  to  create  a  single  transformation  matrix. 

[Transform]  =  [t]  [Pe]  [Pa]  [P-e]  [T-]  (4.2) 


4.6  Hub  Image  Placement 

The  method  of  images  is  used  for  the  hub  and  duct  image  placement  -  as  is  done 
for  the  blade  vortex  lattice  hub  and  duct.  Because  there  is  no  general  guarantee 
that  a  shaft  extends  downstream  of  the  propeller  through  the  length  of  the  wake,  the 
hub  images  are  more  an  accounting  method  to  convect  the  hub  image  vortices  which 
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Point  to  be  imaged 


Point  to  be  imaged 

r_one 

Hubmost  Point 

1  r_zero 

rjmage  ^  |  , 

Generator  Axis 

Profile  View 

Figure  4-5:  Terminology  used  in  determining  the  position  of  the  wake  hub  image 
lattice  point  positions. 


are  required  for  enforcement  of  the  kinematic  boundary  condition  on  the  hub  in  the 
presence  of  the  blade  vortices. 

Due  to  the  non-axisymmetric  nature  of  the  wake  vortex  loop  position,  it  is  not 
possible  to  use  a  common  axis  center  to  place  the  image  vortices.  Rather,  the  position 
of  the  wake  hub  images  must  follow  from  the  position  of  the  wake  vortex  lattice 
elements. 

As  normal,  the  radius  of  the  hub  image  point  is  then  given  by  Equation  4.3. 


image  — 


(rp  -  HGAPf 

n 


(4.3) 


Here  rp  is  the  distance  from  hubmost  wake  lattice  point  to  the  generator  axis  and  ri 
is  the  distance  from  the  point  to  be  imaged  to  generator  axis.  The  tricky  part  is  to 
now  define  the  center  axis  from  which  rp  and  ri  are  measured.  As  luck  would  have 
it,  though,  the  center  axis  was  previously  defined  in  convecting  the  wake  in  the  first 
place. 

The  hub  image  location,  as  shown  in  Figure  4-5  is  placed  in  a  position  coplanar 
with  the  innermost  wake  vortex  lattice  element  to  prevent  undesirable  lattice  skew. 
The  axial  location  of  the  image  vortice  is  set  to  match  the  axial  location  of  the 
hubmost  lattice. 
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Figure  4-6:  Velocity  components  of  a  10  degree  inclined  flow  as  seen  by  a  propeller. 
Plots  of  the  axial,  radial,  and  tangential  component  of  the  inflow  at  the  propeller 
seven-tenths  radius. 


Force  Component 

Aligned  Wake 

Non-aligned  Wake 

A% 

Steady  Kt 

0.15799 

0.15740 

-0.37 

Steady  IOA'q 

0.29348 

0.29443 

0.32 

1st  Harmonic  Axial  Force 

0.01849 

0.01458 

21 

1st  Harmonic  Vertical  Force 

0.00158 

0.00217 

37 

1st  Harmonic  Side  Force 

0.01055 

0.00823 

22 

Table  4.1:  Aligned  and  Non-aligned  wake  results  for  10  degree  angle  of  inclination. 

4.7  Aligned  Wakes  with  Prescribed  Effective  In¬ 
flow 

The  results  presented  here  represent  calculations  made  with  a  prescribed  effective 
inflow  inclined  at  10  degrees.  Fully  coupled  results  are  presented  in  Chapter  6.  A 
standard  three-bladed  propeller,  DTMB  Propeller  4119,  is  used  for  all  calculations. 

The  first  result  is  to  show  the  wake  position  when  aligned  with  a  ten  degree 
angled  inflow  as  shown  in  Figure  4-2.  The  axial,radial,  and  tangential  components  of 
the  inflow  are  shown  in  Figure  4-6.  The  once  per  revolution  variation  in  tangential 
velocity  component  will  produce  a  wake  with  non-constant  pitch  angle.  The  resulting 
differences  in  predicted  operating  conditions  between  using  the  aligned  wake  model 
and  the  non-aligned  wake  model  are  shown  in  Table  4.1.  The  differences  in  lateral 
forces  are  maily  due  to  phase  differences.  Notice  that  while  the  mean  forces  are 
generally  similar,  the  predicted  first  harmonic  forces  differ  by  an  order  of  magnitude 
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Circulation  at  the  0.7  Radius  in  10  Degree  Inclined  Inflow 


Figure  4-7:  With  the  propeller  operating  in  a  10  degree  inclined  inflow,  this  plot 
shows  the  circulation  at  the  ^  radius  as  the  propeller  rotates  through  360  degrees. 

Figure  4-7  highlights  the  difference  in  resulting  blade  solution  for  aligned  and 
non-aligned  wakes.  Figure  4-7  shows  the  blade  circulation  at  the  radius  as  the 
propeller  sweeps  through  360  degrees  in  a  10  degree  inclined  flow.  Note  that  the  mean, 
or  steady  circulation,  is  nearly  the  same,  but  the  prediction  in  the  V*  harmonic  shows 
phase  and  amplitude  differences. 


4.7.1  Once  per  Revolution  Change  in  Axial  Velocity 

The  second  test  is  to  vary  the  wake  pitch  angle  by  creating  a  stylized  flowfield  with 
varying  axial  velocity.  The  components  of  the  flowfield  are  shown  in  Figure  4-8.  This 
is  a  simplified  model  for  the  effects  of  upstream  appendages.  The  variation  in  axial 
velocity  will  produce  a  varying  wake  pitch  angle.  Again,  the  propeller  program  is  run 
with  and  without  wake  alignment. 

The  results  of  this  analysis  are  shown  in  Figure  4-9.  Because  there  is  negligible 
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Force  Component 

Aligned  Wake 

Non-aligned  Wake 

A% 

Steady  Fx 

-0.04874 

-0.05068 

4.00 

Steady  Fy 

0.00274 

0.00270 

1.50 

Steady  Fz 

0.02822 

0.02891 

2.50 

1st  Harmonic  Fx 

0.01555 

0.01311 

15.7 

1st  Harmonic  Fy 

0.00214 

0.00296 

38.3 

1st  Harmonic  Fz 

0.00752 

0.00574 

23.7 

Table  4.2:  Unsteady  blade  forces  produced  by  a  propeller  in  an  axially  non-uniform 
flow  with  and  without  wake  alignment.  There  are  large  penalties  associated  with  not 
aligning  the  wake. 


Axial  Velocity  ^  Radial  Velocity  ,,  Tangential  Velocity 


Figure  4-8:  Velocity  components  for  a  once  per  revolution  variation  in  the  magni¬ 
tude  of  the  axial  component  of  the  inflow.  Plots  of  the  axial,  radial,  and  tangential 
component  of  the  inflow  at  the  propeller  seven- tenths  radius. 


steady  offset,  the  mean  predicted  propeller  Kt  and  IOKq  differ  insignificantly.  The 
once  per  revolution  variation,  though,  highlights  the  difference  in  unsteady  perfor¬ 
mance.  Table  4.7.1  shows  the  differences  in  1st  harmonic  forces. 


4.8  Wake  Lattice  Spacing 

Due  to  the  requirement  that  the  vorticity  present  on  a  given  wake  panel  convects 
downstream  in  the  time  interval  between  successive  blade  clicks,  the  size  of  the  wake 
vortex  lattice  loops  are  directly  determined  by  this  time  step  (Number  of  Time  Steps 
per  revolution  in  PUF  terminology).  Because  each  shed  vortex  loop  corresponds 
to  the  distance  between  angular  blade  positions,  a  large  number  of  blade  steps  per 
revolution  is  required  to  keep  a  sufficiently  small  wake  discretization.  But  an  increased 
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Circulation  at  the  7/1 0th  radius 


Figure  4-9:  4119  circulation  at  the  seven-tenths  radius  for  axially  non-uniform  inflow. 


number  of  time  steps  requires  longer  solution  times.  This  problem  is  highlighted  by 
the  required  angular  resolution  which  follows  directly  from  the  user-specified  number 
of  time  steps  per  revolution  (NTSR). 


The  propeller  time  step  is  given  by: 


^Tpj-0p  — 


2J 

NTSR 


(4.4) 


Which  directly  leads  to  the  following  angular  increment  between  each  successive  blade 
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(4.5) 


Notice  that  the  angular  discretization  is  independent  of  advance  coefficient,  J,  and 
depends  only  upon  the  user  input  number  of  blade  positions  per  revolution.  This 
means  that  if  it  is  possible  to  create  a  wake  discretization  scheme  independent  of 
blade  time  steps  it  will  also  be  problem  independent.  This  would  obviously  be  a  more 
robust  algorithm  in  general. 

This  discretization  error  is  seen  in  the  convergence  rate  of  the  blade  forces  with 
changes  in  pseudo  time  step  (the  NTSR  variable  using  PUF  terminology)  shown  in 
Figure  4-10.  This  test  was  conducted  using  a  uniform  effective  inflow  and  maintaining 
constant  lattice  size  on  the  blade.  These  results  are  similar  to  those  presented  by 
Warren  [30].  Under  the  older  wake  growing  scheme,  where  wake  loop  size  is  dependent 
upon  time  step  discretization,  the  convergence  of  Kt  and  Kq  such  that  they  are  no 
longer  dependent  upon  NTSR  comes  with  a  stiff  penalty  in  the  form  of  computational 
expense.  Notice  from  Figure  4-10  that  the  wake  growing  scheme  is  unstable  at  low 
NTSR.  A  goal  of  this  research  is  to  remove  the  solution  dependence  on  NTSR  while 
maintaining  low  computation  times. 


4.9  Higher  Order  Wakes 

There  are  two  methods  which  could  be  used  to  remove  the  dependence  between  the 
propeller  discrete  time  stepping  and  the  wake  loop  discretization. 
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Figure  4-10:  The  convergence  of  the  propeller  solution  is  directly  determined  by  the 
pseudo  time  stepping,  which  is  controlled  via  the  user-defined  Number  of  Time  Steps 
per  Revolution  (NTSR). 

1.  Higher  order  vortex  strength  representation.  The  coarse  discretization  could  be 
overcome  with  higher  order  representation  of  the  vortex  loop  strengths. 

2.  Higher  order  panels.  If  the  shed  vortex  loops  can  be  considered  constant 
strength  singularity  panels,  then  a  higher  order  panel  would  accurately  repre¬ 
sent  the  geometry,  while  minimizing  the  computational  expense.  The  curvature 
effects  are  missed  for  very  coarse  panels. 

These  methods  directly  result  from  the  fact  the  influence  velocity  of  the  wake  is  due 
to  the  product  of  the  influence  coefficient  and  the  vortex  loop  strength 

=  ^n[WIF]i,k  (4.6) 

k=l 

The  secondr  method  is  the  better  choice,  since  a  higher  order  panel  will  also  solve 
the  instability  in  the  vortex  loop  placement  scheme  associated  with  the  current  lower 
order  scheme. 
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4.9.1  Higher  Order  Wake  Loop  Scheme 


A  classical  higher  order  loop  would  be  created  by  defining  the  perimeter  of  the  loop, 
and  using  a  higher-order  influence  function  to  return  the  influence  of  that  loop.  The 
degree  would  depend,  then,  on  the  higher  order  derivatives  at  the  wake  loop  edges. 
This  is  nearly  equivalent  to  breaking  the  loop  into  smaller  loops,  calculating  the 
influence  of  each  smaller  loop,  and  re-assembling  these  smaller  influences  into  an 
overall  influence  coefficient  for  the  total  loop.  Recall  that  global  loop  size  is  still 
dictated  by  the  propeller  discrete  time  stepping  (since  shed  vorticity  from  the  blades 
at  each  pseudo  time  step  must  be  convected  in  the  wake).  The  advantages  of  the  loop 
subdivision  scheme  are, 

•  wake  curvature  effects  are  captured 

•  more  stable  than  large  stepping  scheme  currently  used 

•  ability  to  easily  integrate  a  streamwise  strength  variation  scheme 

Figure  4-11  shows  pictorially  the  resulting  wake  loop  discretization  with  and  with¬ 
out  subidivisions.  The  smaller  subdivided  panel  better  capture  the  real  curvature 
effects  of  the  wake,  and  reduce  the  mis-alignment  between  the  wake  panel  tengency 
and  the  local  fluid  flow. 

4.9.2  Convergence  with  Subdivided  Wake  Loops 

Figure  4-12  shows  the  results  for  varying  the  subdivisions  of  each  wake  vortex  loop. 
Notice  that  KT  and  lOKQ  converge  after  the  wake  loops  are  subdivided  to  an  equiva¬ 
lent  NTSR=60  level,  and  requires  one-third  of  the  computation  time.  Based  on  these 
results,  the  algorithm  is  set  to  subdivide  the  wake  vortex  loops  to  use  an  equivalent 
pseudo  timestepping  which  would  result  if  the  user  had  specified  NT  SR  =  60  in 
the  old  scheme.  The  pseudo  time  stepping  in  the  wake,  then,  is  not  equivalent  to 
Equation  4.4  but  is  given  by  Equation  4.7. 
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Figure  4-11:  Subdividing  the  loop  sizes  required  by  the  pseudo  time  step  into  smaller 
loops  creates  a  higher  order  influence  panel.  The  solid  dashed  lines  are  the  actual 
wake  loop  sizes  as  dictated  by  the  pseudo  time  step  and  the  lighter  solid  lines  are  the 
subdivided  wake  loop  locations.  There  are  five  subdivisions  per  time  step. 


(4.7) 

And  the  equivalent  rotational  component  added  at  each  time  step  is  shown  in  Equa¬ 
tion  4.8. 


^Prop  =  ^  (4.8) 

The  results  show  a  70  percent  reduction  in  computation  time  for  equivalent  accu¬ 
racy  and  convergence.  Computation  time  could  be  further  reduced  by  implementing 
a  higher  order  influence  calculation  since  the  computation  burden  is  now  evenly  di¬ 
vided  between  the  wake  growing  algorithm  using  small  pseudo  time  steps  and  the 
calculation  of  the  influences  for  each  of  these  smaller  loops.  The  time  savings  results 
from  many  factors: 


57 


•  reduced  number  of  angular  positions  requires  fewer  time  steps  to  convect  the 
starting  vortex  downstream 

•  the  force  calculation  directory  depends  on  the  number  of  blade  panel  positions 

•  less  memory  required 

The  use  of  the  user-defined  parameter  for  the  number  of  timesteps  per  revolution 
should  now  only  be  set  higher  than  the  number  of  blades  to  capture  unsteady,  higher 
harmonic  effects.  There  are  many  calculations  which  do  not  require  the  calculation 
of  higher  harmonic  unsteady  forces,  which  will  benefit  from  the  70  percent  speed-up. 

•  Maneuvering  vehicle  calculations  typically  only  require  1st  and  2nd  harmonic 
unsteady  forces. 

•  When  coupled  with  an  external  flow  solver,  such  as  a  RANS  code,  it  does  make 
any  sense  at  all  to  try  and  calculate  harmonic  forces  higher  than  the  azimuthal 
grid  resolution.  For  instance,  given  a  seven  bladed  propeller  coupled  with  a  , 
RANS  code  with  36  grid  cells  in  the  azimuthal  direction,  the  user  should  not 
expect  to  resolve  any  harmonics  higher  than  2. 

4.9.3  Solution  Validation  with  Subdivided  Wake  Loops 

The  numerical  accuracy  of  the  new  wake  subdivision  scheme  is  found  by  comparing 
the  results  of  PUF  with  and  without  the  wake  subdivision  scheme,  where  there  are  a 
large  number  of  wake  loops  in  the  non-subdivided  scheme.  This  check  is  necasarray 
due  to  the  large  scale  changes  in  array  indexing,  solution  advancement,  and  influence 
looping  the  wake  subdivision  scheme  required. 
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Figure  4-12:  Using  a  pseudo  higher  order  influence  scheme  through  subdividing  the 
wake  loops  results  in  less  computation  effort  for  convergence. 

The  numerical  results  for  circulation  distribution  with  the  wake  subdivision  scheme 
is  shown  in  Figure  4-13.  The  integrated  results  for  Kt  and  \^Kq  are, 

NTSR  Kt  \0Kq 

PUF14^  with  subdivided  wakes  3  0.15173  0.28532 

PUF14^  without  subdivided  wakes  60  0.15174  0.28534 

PUF14^  without  subdivided  wakes  60  0.15231  0.28184 

^  Individual  Wakes 
^  Imaged  Wakes 

It  is  hypothesized  that  this  difference  is  due  to  round  off  errors  present  when  rotating 
the  imaged  wakes  to  the  correct  angular  positions. 

4.10  Validations  of  Aligned,  Subdivided  Wake  Scheme 

To  test  the  subdivision  scheme,  Propeller  4119  is  given  an  unsteady  effective  inflow. 
While  the  blade  lattice  size  is  held  flxed  at  15X15,  the  number  of  discrete  time  steps 
is  increased  to  the  maximum  memory  of  the  machine.  The  results  are  presented 
in  Table  4.3.  Notice  that  the  results  do  vary  slightly.  This  is  due  to  the  different 
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Figure  4-13:  Validation  of  PUF14  using  the  wake  subdivision  scheme  with  PUF14 
without  the  wake  subdivision  scheme  in  a  uniform,  axial  inflow. 

NTSR  Kt  IOKq 

15  0.15380  0.28503 

30  0.15423  0.28643 

45  0.15453  0.28731 

Table  4.3:  Affect  of  NTSR  on  subdivided  wake  scheme  results.  The  difference  of  less 
than  1%  is  attributable  to  the  variable  spacing  in  the  first  wake  panel,  which  does 
depend  upon  NTSR 


subdivision  sizes  present  in  the  first  wake  lattice. 


4.10.1  Convergence  Test  of  Aligned,  Subdivided  Wake  Scheme 

The  first  convergence  test  is  to  investigate  if  the  propeller  solution  with  the  subdivided 
wakes  requires  a  finer  blade  grid  than  the  constant-sized  wake  panel  solution.  To  make 
the  investigation  more  interesting,  the  propeller  code  with  subdivided  wakes  is  tested 
against  the  propeller  code  with  constant-sized  wakes  in  a  non-uniform  inflow  which 
is  specified  as  the  effective  flowfield.  The  convergence  rate  with  blade  lattice  size  is 
shown  in  Figure  4-14.  These  results  indicate  that  the  users  of  the  new  code  do  not 
need  to  change  their  thought  process  on  the  correct  settings  for  blade  discretization. 
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Discretization  Convergence  Tests 
Subdivided  and  Non-subdivded  Wake  Loops 


Figure  4-14:  Convergence  rate  of  the  new  wake  discretization  scheme  in  a  non-uniform 
inflow.  Compared  to  the  previous  scheme  of  large,  non-subdivided  wake  panels,  the 
convergence  rate  is  similar  with  changes  in  blade  lattice  density. 


61 


0.34 

0.32 

0.3 

0.28 

0.26 

0.24 

0.22 

0.2 

0.18 

0.16 

0.14 

0.12 

0.1 


Discretization  Convergence  Tests  in  PUF  with 
Fully  Aligned,  Subdivided  Wakes  in  Uniform 
and  Non-uniform  Inflows 


Propeller  4119 

Uniform  Inflow  Vel,^  =  1.0 

Non  Uniform  Inflow  with  VeL  =  l+sin(2  theta) 
IHUB=1/2NKEY 

-  ^ 

— V - V - w - V 

r  i  1  »  -■  r  I  I  t _ I _ |_J _ I _ I  I  I  I _ i _ _ I _ I _ ] _ I _ |__| _ I 

5  10  15  20  25  30 

Blade  Lattice  Size  (NxM) 


Figure  4-15:  Mean  Kt  and  IOjFsTq  results  are  shown  for  PUF  using  fully  aligned, 
subdivided  wakes  in  both  uniform  and  non-uniform  inflows.  Grid  convergence  occurs 
at  a  lattice  size  of  around  20. 

Mean  Kt  and  lOi^g  results  are  shown  in  Figure  4-15  for  PUF  using  fully  aligned, 
subdivided  wakes  in  both  uniform  and  non-uniform  inflows.  These  results  indicate 
that  grid  convergence  occurs  at  a  lattice  size  of  around  20. 


4.10.2  Propeller  4119  in  Non  Uniform  Axial  Inflow 

The  varying  pitch  angle  of  the  wake  as  shown  in  Figure  4-1  is  compared  with  the 
wake  grown  in  previous  methods  in  Figure  4-16.  While  there  is  some  difference  in  the 
performance  prediction  between  the  aligned  and  non-aligned  calculations,  the  largest 
difference  shows  up  in  the  prediction  of  the  blade  harmonic  forces. 
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Figure  4-16:  A  non-uniform  inflow  directly  affects  the  pitch  angle  of  the  wake  sheets. 
Here  the  non-uniformity  is  produced  by  a  step  function  change  in  the  axial  velocity. 
Aligning  the  wakes  with  the  local  flow  produces  radical  changes  in  the  wake  pitch 
angle.  The  aligned  wake  is  the  same  as  shown  in  Figure  4-1. 
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Chapter  5 

Coupled  Algorithms 


The  major  shortcoming  of  the  lifting  surface  formulation  described  above  is  the  in¬ 
ability  to  model  the  vorticity  interaction.  As  previously  mentioned,  the  propeller 
induction  velocities  cause  a  redistribution  of  the  vorticity  present  in  the  inflow.  This 
interaction  is  critical  in  the  design  of  modern,  full  stern  submarines,  where  the  stern 
prismatic  coefficient  is  so  high  that  the  propulsor  actually  inhibits  separation  of  the 
boundary  layer. 

The  second  area  of  propulsor  design  affected  by  this  shortcoming  is  in  the  design 
of  multi-element  propulsors.  The  problem  is  that  at  points  where  the  trailing  wake 
vorticity  from  the  upstream  blade  rows  intersect  the  control  points  of  the  downstream 
blades,  the  mathematics  of  the  problem  produces  a  singular  solution.  Multi-element 
propulsor  design  is  important  in  the  design  of  modern  naval  combatant  propulsors, 
and  in  the  design  of  water] et  pumps. 


5.1  Time- averaged  Induced  Velocity 

When  coupling  a  potential-based  lifting-surface  blade  solver  with  a  RANS  code,  the 
need  arises  to  correctly  evaluate  the  effective  velocity  given  the  total  velocity  solved 
for  by  the  RANS  code. 
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In  the  propeller  code,  the  velocity  at  any  field  point  can  be  broken  down  as  follow: 


V®  =  v«,  +  V,®  +  K® 


(5.1) 


where:  is  the  effective  inflow 

is  the  time-averaged  induced  velocity 
y®  is  the  fluctuating,  local  blade-to-blade  induced  velocity 
Since  the  effective  inflow  is  not  measurable  in  the  RANS  code,  the  assumption  is 
made  that  the  mean  of  the  local  blade-to-blade  flow  is  zero, 

y®  =  0  (5.2) 

This  then  directly  leads  to  the  final  statement  that  the  time-averaged  effective  velocity 
is  equivalent  in  the  RANS  code  ,  as  indicated  by  the  y®  notation  and  the  propeller 
code,  such  that  the  effective  velocity  is  just  the  total  velocity  minus  the  time-averaged 
induced  velocity. 


K?/  =  y?,/ = 


(6.3) 


This  statement  ignores  the  effect  of  the  fluctuating  blade-to-blade  component  on  the 
vorticity  redistribution,  explained  in  Section  B.l,  but  is  countered  by  the  physical 
observation  that  the  higher  harmonically  fluctuating  blade  to  blade  flows  quickly 
decay  away  from  the  blades. 

The  computational  burden  then  comes  from  the  requirement  to  compute  the  time- 
averaged  induced  velocity  at  all  boundary  condition  points.  In  its  most  elemental 
form,  this  calculation  involves  determining  the  velocity  induced  at  a  field  point  in 
space  by  a  large  set  of  three  dimensional  vortex  segments  rotated  about  an  axis  to 
form  conic  vortex  sections,  where  the  strength  of  the  vorticity  varies  with  respect  to 
angular  position  on  the  conic. 
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5.2  Velocity  Induced  by  a  Vortex  Segment 


The  velocity  induced  by  a  vortex  segment  on  a  point  in  space  is  known  precisely  from 
the  Biot-Savart  Law  [14], 


dv  = 


r  Rx  ds 
Att 


(5.4) 


where:  dv  is  the  induced  velocity  vector 

R  is  the  vector  from  the  vortex  element  to  the  field  point 
ds  is  the  unit  vector  along  the  vortex  segment 
R^  is  the  cube  of  the  distance  vector 

Integrating  Equation  5.4  along  the  length  of  the  segment  gives  the  total  induced  veloc¬ 
ity  upon  the  field  point.  The  numerical  implementation  of  Equation  5.4  is  described 
more  fully  by  Greeley  [7]. 


5.2.1  Discrete  Calculation  of  the  Velocity  Induced  by  a  Conic 
Section 


Using  the  Biot-Savart  law  as  implemented  numerically  by  Greeley  [7]  naturally  leads 
to  representing  the  conic  section  as  a  series  of  vortex  sticks  swung  about  the  rc-axis. 
This  is  shown  in  Figure  5-1.  Because  the  strength  of  the  vortex  segment  varies  with 
angular  position,  it  follows  that  the  induced  velocity  of  the  conic  section  can  be  found 
through  integrating  the  effect  of  the  discrete  vortex  segments. 


V  = 


^  r(^)  Rxds 
47r  R^ 


de 


(5.5) 
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Single  Vortex  Segment  in  Free  Space . rotated  about  the  x-axis  at  discrete  angular  positions 


Figure  5-1:  Conic  section  made  up  of  discrete  vortex  segments.  This  illustration 
shows  how  discrete  vortex  segments  can  be  combined  to  form  a  nearly  continuous 
sheet. 


5.3  Velocity  Induced  by  a  Ring  Vortex  of  Harmon¬ 
ically  Varying  Strength 


An  alternative  method,  presented  by  Buchoux  [1]  for  the  uniform  axisymmetric  case, 
and  refined  by  Greeley  is  to  represent  the  conic  section  as  rings  of  vorticity.  The  Biot- 
Savart  equation  is  still  the  governing  equation,  but  instead  of  integrating  around 
0  to  27r,  vortex  rings  are  integrated  along  the  surface  of  the  conic  section.  The 
derivation  for  the  induced  velocity  due  to  a  vortex  ring  of  uniform  strength  is  shown 
in  Kucheman  [18].  The  final  resulting  equation  for  the  velocity  induced  by  a  vortex 
ring  of  harmonically  varying  strength  was  shown  first  by  Greeley  [6]  who  implemented 
several  computer  codes  to  solve  the  model  problem. 


5.3.1  Velocity  Induced  by  Constant  Strength  Vortex  Ring 


As  shown  in  Kucheman  [18],  a  constant  strength  vortex  ring  induces  only  an  axial 
and  radial  component  of  velocity.  The  x  component  of  velocity  is  given  by  integrating 
the  ring  around  the  circumference, 


Vx{x,  r) 


r  rcos((^  -(j)')- I 

47rr'  Jq  +  r^  +  1  —  2rcos{4>  —  (t>'Y 


(5.6) 
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This  integral  is  expressed  in  terms  of  elliptic  integrals  through, 


Vx{x,r)  = 


47rr'  yjx^  {r  +  +  (r  —  1)^] 


h 


where 


h 


p2'K 

Jo 


rcos{4>  —  (/>')  —  1 


_  2r[co5 (0-00+1] 


x2+(r+l)2 


1  + 


2r 


1  — cos(0— 0')j 


®2+(r— 1)2 


/i  can  be  further  simplified  by  letting, 


x^  +  {r  +  1)2 

such  that. 

This  gives  the  final  result  that, 

r  1 

=  - - ^===K(k)  - 

2nr  +  (r  +  1)2 

A  similar  derivation  can  be  followed  to  arrive  at  the  radial  velocity  induced  by  a 
vortex  ring.  Due  to  the  axial  symmetry  of  the  problem  there  can  be  no  tangential 
component  of  the  induced  velocity 


1  + 


2(r  -  1) 
a:2  +  (r  —  1)2  J 


E{k) 


(5.7) 


5.3.2  Non-constant  Strength  Vortex  Ring  Induced  Velocity 

Non-constant  strength  vortex  rings  can  be  modeled  through  replacing  T  in  Equation 
5.6  with  r„sm(n^).  Modeling  the  various  harmonic  component  of  T,  Equation  5.6  is 
modified  to, 

V.(x,r)  =  -^  r 

47rr'  Jq  Y^a;2  -1-  r2  +  1  —  2rcos{(j)  —  4>')^ 

Equation  5.8  can  again  be  reduced  to  a  complete  elliptical  integral  of  the  third  kind. 
All  three  components  of  velocity  are  present,  now,  due  to  the  loss  of  axial  symmetry. 
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5.3.3  Conic  Sections  by  Successively  Integrating  Vortex  Rings 

The  conic  section  problem  is  then  solved  by  successively  integrating  the  vortex  rings 
along  the  length  of  the  conic,  and  properly  accounting  for  the  phase  shift  and  ampli¬ 
tudes  of  the  induced  velocities.  Due  to  the  linear  nature  of  the  problem,  the  solution 
is  sped  up  by  finding  the  complex  response  phase  and  amplitude  response  at  a  point 
defined  by  its  axial  and  radial  coordinates.  The  response  at  other  points  along  the 
same  radial  ring  are  phase  shifted  from  this  initial  solution. 

Adaptive  Quadrature  Method 

An  adaptive  quadrature  scheme  is  required  due  to  the  highly  non-linear  nature  of 
the  influence  velocity  calculation.  The  Romberg  integration  scheme  [26]  with  the 
Richardson  extrapolation  method  [26]  was  implemented.  The  basic  methodology  is 
outlined  in  Appendix  C.l. 

At  each  subdivision  step,  a  trapezoidal  integration  is  performed  followed  by  the 
Richardson  error  extrapolation.  Based  on  the  convergence  check,  a  further  subdivision 
is  performed  unless  the  maximum  subdivision  level  has  been  reached. 

5.3.4  Discrete  and  Analytic  Model  Comparison 

The  necessity  of  the  analytic  model  can  be  shown  through  comparing  the  discrete 
and  analytic  induced  velocities  calculated  for  simple  model  lamp  shade  problems. 
The  lamp  shade  is  formed  by  rotating  a  three  dimensional  element  of  vorticity  about 
the  a;-axis  as  seen  in  Figure  5-1.  To  simulate  the  harmonic  effects  found  in  unsteady 
propeller  forces,  harmonically  varying  strengths  are  assigned  to  the  elemental  vortex 
as  it  traverses  360  degrees,  where  the  strength  of  the  segment  is  dependent  upon  the 
angular  position  of  the  element. 

The  question  as  to  how  fine  a  discretization  level  is  required  to  produce  results 
comparable  to  the  analytic  results  must  be  answered.  This  is  done  through  two  model 
problems  where  a  field  point  is  moved  relative  to  a  conic  section  represented  by  sticks 
of  vortex  segments. 
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Figure  5-2:  Two  model  problems  to  highlight  the  discretization  required  to  reproduce 
analytic  results. 


The  first  model  problem  involves  forming  a  conic  section  of  straight  vortex  seg¬ 
ments  and  rotating  a  field  point  around  the  outside  of  the  conic  section.  Increasing 
the  number  of  discrete  segments  makes  the  resulting  induced  velocity  calculations 
converge  to  the  analytic  results  as  shown  in  Figure  5-3.  The  method  of  using  dis¬ 
cretized  lateral  elements  does  not  converge  until  a  high  number  of  elements  are  used, 
while  the  method  of  successively  integrating  analytic  rings  converges  quite  quickly. 
Note  the  instability  in  the  predicted  induced  velocities  for  a  medium  resolution  level 
of  discretization. 

The  second  validation  is  to  again  calculate  the  influence  of  a  conical  section  of  vor- 
ticity  upon  a  field  point  as  the  field  point  approaches,  and  passes  through,  the  lamp 
shade.  Again,  the  conical  section  has  non-uniform,  harmonically  varying  strength. 
The  results  shown  in  Figure  5-4  show  that  a  high  number  of  discretizations  are  needed 
for  the  discrete  solution  to  approach  the  analytic  solution.  The  method  of  using  dis¬ 
cretized  lateral  elements  does  not  converge  until  a  high  number  of  elements  are  used, 
while  the  method  of  successively  integrating  analytic  rings  converges  quite  quickly. 
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Figure  5-3:  Induced  velocity  upon  a  field  point  as  it  traverses  around  a  conical  vortex 
sheet  as  predicted  by  both  the  analytic  method  and  discretized  stick  elements. 


72 


Distance  of  Field  Point  from  Lampshade  Surface 


Figure  5-4:  Induced  velocity  upon  a  field  point  as  it  pierces  a  conical  vortex  sheet  as 
predicted  by  both  the  analytic  method  and  discretized  stick  elements. 


5.4  Induced  Velocity  Methodology  Implementation 


The  analytic  expression  for  the  velocity  induced  by  a  non-uniform  vortex  ring  was 
implemented  as  part  of  a  new  fully  unsteady  three  dimensional  vortex-lattice  lifting- 
surface  code,  PUF-15.  The  basic  procedure  follows  as, 

1.  At  the  end  of  the  blade  solution,  calculate  blade  circulation  harmonic  amplitudes 
and  phases. 

2.  For  each  blade  loop,  calculate  the  analytic  induced  velocity  based  on  an  adap¬ 
tive  Romberg  iteration  scheme  which  calls  the  harmonic  vortex  ring  influence 
routines. 

3.  The  wakes  are  handled  as  axisymmetric  annuluses,  similar  to  the  blades,  with 
the  proper  shed  vortex  strength. 

4.  The  induced  velocity  is  then  the  recombination  of  the  influence  and  the  circu¬ 
lation  amplitude,  offset  by  the  phase  difference  between  the  response  and  the 
input. 

5.4.1  Numerical  Validations 

Axisymmetric  ring  vortices  of  constant  strength  were  shown  by  Buechoux  [1]  to  have 
an  analytic  solution,  and  an  induced  velocity  routine  based  on  those  derivations  was 
implemented  by  Bechoux  in  the  steady  propeller  analysis  code,  PBD-14.  By  ana¬ 
lyzing  the  unsteady  propeller  in  a  steady  flowfleld  the  results  from  the  two  analysis 
techniques  are  demonstrable. 

Figure  5-5  shows  the  comparison  of  the  analytic  routines  implemented  by  Be¬ 
choux  in  the  steady  propeller  code  PBD  and  the  analytic  routines  implemented  in 
the  unsteady  propeller  program  PUF15  by  this  research.  Both  codes  use  an  analyt¬ 
ical  representation  of  the  vorticity  as  vortex  rings  and  solve  the  problem  via  elliptic 
integrals.  The  good  agreement  and  solution  smoothness  over  the  entire  blade  val¬ 
idates  the  new  approach.  Comparison  between  the  two  analytic  codes  shows  very 
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Figure  5-5:  Comparison  of  induced  velocity  calculated  in  steady  propeller  analysis 
program,  solid  contour  lines,  and  the  unsteady  propeller  analysis  program  PUF  in 
dashed  lines. 


good  agreement.  Within  an  unsteady  code,  a  comparison  can  be  made  between  the 
induced  velocities  produced  by  the  discrete  algorithm  and  the  new  analytic  algorithm. 
Comparison  between  the  discrete  modeling  algorithm  implemented  in  PUF14  and  the 
analytic  algorithm  implemented  in  PUF15  is  shown  in  Figure  5-6.  The  largest  dis¬ 
crepancies  are  found  near  the  tip,  where  the  blade  vortex  lattice  segments  are  densely 
spaced.  This  discrepancy, however,  is  balanced  by  the  fact  that  there  is  very  little 
load  at  the  tip,  which  minimizes  the  error  introduced  into  the  coupling. 


5.5  Grid  Intersection  Methodology 

The  second  major  piece  in  coupling  the  vortex-lattice  lifting-surface  code  and  a  RANS 
code  is  to  interpolate  the  velocities  from  the  RANS  grid  to  the  propeller  grid,  and  to 
interpolate  the  propeller-generated  forces  within  the  swept  volume  of  the  propeller  to 
the  proper  RANS  cells. 

The  interpolation  of  continuous  propeller  forces  to  the  RANS  grid  is  the  more 
complex  of  the  two  interpolations  since  volume  interpolations,  vice  point  interpola¬ 
tions,  must  be  done.  Previous  coupling  methods  have  used  an  exhaustive  polyhedron 
overlap  search  algorithm  which  finds  the  percent  overlap  between  each  PUF  cell  and 
each  RANS  grid  cell. 
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Figure  5-6:  Comparison  of  induced  velocity  calculated  in  unsteady  propeller  analysis 
code  via  discrete  (solid  lines)  and  analytic  (dashed  lines)  methods.  The  comparison 
of  induced  velocity  vectors  in  (d)  shows  the  large  discrepancy  at  the  tip. 
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Figure  5-7:  Grid  subdivision  to  transfer  forces  between  PUF  and  RANS.  The  upper 
figure  is  the  circumferential  PUF  blade  grid  where  body  forces  are  know.  The  bottom 
grid  is  the  actual  RANS  grid,  and  the  middle  grid  is  the  subdivided  RANS  grid,  which 
is  the  fine  transfer  mesh. 

This  thesis  has  implemented  a  new  algorithm  which  uses  a  cell  subdivision  tech¬ 
nique.  The  process  is  to  subdivide  the  RANS  cells  into  smaller  cells.  It  is  then  a 
very  fast  search  to  determine  if  the  center  of  the  subdivided  RANS  cell  lies  within 
any  of  the  PUF  cells.  If  the  center  of  the  subdivided  RANS  cell  is  found  to  lie  within 
the  swept  volume  of  the  propeller,  a  linear  interpolation  is  done  from  the  edges  of 
the  PUF  cell  to  place  the  correct  force  density  in  the  RANS  cell.  Recall  that  RANS 
solvers  introduce  the  propeller  forces  as  source  terms  on  the  right  hand  side  of  the 
RANS  equation  [29].  This  formulation  also  leads  to  the  propeller  forces  being  called 
body  forces.  A  pictoral  of  this  process  is  shown  in  Figure  5-7. 

Numerical  experiments  show  that  subdividing  the  RANS  cells  by  a  factor  of  10 
along  each  side,  such  that  there  are  1, 000  subdivided  cells  within  each  RANS  cell, 
keeps  the  interpolation  error  well  under  0.10%.  And  because  the  searching  algorithm 
now  involves  the  search  for  a  point  within  a  limited  set  of  domains,  vice  polyhedron 
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overlapping  volumes,  the  algorithm  is  actually  faster. 


i 
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Chapter  6 


Coupled  Validations 


Coupling  the  propeller  and  RANS  codes  allows  for  numerical  validation,  as  well  as 
validation  against  experimental  test  results. 

Section  6.1  presents  results  for  coupling  the  propeller  code  PUF-15  with  two  sep¬ 
arate  RANS  codes,  IFLOW-41  and  UNCLE-3D.  Two  RANS  codes  are  used  to  assess 
the  relative  errors  due  to  the  discretizations  and  formulations  of  the  propeller  code, 
the  coupling  code,  and  the  RANS  code.  Section  6.2  presents  the  results  for  the  same 
propeller  4119  inclined  at  10  degrees  with  respect  to  the  horizontal.  Section  6.3 
presents  coupled  results  for  the  DTMB  Propeller  4679  driven  by  a  downstream  shaft 
inclined  at  7.5  degrees.  Experimental  data  recorded  for  propeller  4679  is  compared 
with  the  predictions  from  this  thesis. 

Figure  6-1  gives  an  overview  of  the  code  interaction,  and  workflow  followed  to 
show  validated  results.  Among  the  four  coupled  cases  presented  here,  all  validations 
were  performed  to  the  greatest  extent  possible.  The  blocks  in  the  center  of  Figure 
6-1  represent  the  codes  used  for  the  analysis,  and  the  dashed  line  encompasses  newly 
created  analysis  codes 

6.1  Propeller  4119  Validations 

The  DTMB  propeller  4119  shown  in  Figure  6-2  represents  nearly  an  ideal  propeller 
to  validate  against  due  to  the  large,  benign  blade  sections,  with  no  rake  and  skew 
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Figure  6-1:  This  schematic  gives  an  overview  of  the  coupled  RANS  cases  run,  and 
the  types  of  analysis  performed.. 


Profile  View 


Figure  6-2;  These  plots  show  the  geometric  simplicity  of  the  DTMB  propeller  4119. 

[4].  Experimental  data  is  available  for  boundary  layer  profiles,  pressure  distributions, 
drag  coefficients,  and  downstream  wake  survey  [12].  This  fact  and  the  geometric 
simplicity  make  4119  an  ideal  validation  candidate  for  numerical  propeller  schemes. 


6.1.1  UNCLE-3D  Coupling  with  Propeller  4119 

The  UNCLE1-3D  RANS  solver  represents  a  modern  day  RANS  solver.  The  computa¬ 
tional  RANS  grid  is  shown  in  Figure  6-3.  It  is  an  0-grid  topology  of  size  75X35X57. 
A  grid  convergence  study  was  not  completed  since  this  validation  work  follows  from 
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Figure  6-3:  Numerical  grid  used  in  the  UNCLE1-3D  RANS  code 


previous  validation  work  [30].  Note  the  great  number  of  cells  across  the  chord  of  the 
propeller.  Experience  has  shown  that  greater  than  10  cells  are  required  across  the 
chord  to  accurate  capture  the  changes  which  occur  in  the  blade  passage. 

Induced  Velocity  Calculations 

Coupling  the  propeller  and  RANS  code  in  a  uniform  inflow,  with  a  no-slip  boundary 
condition  on  the  shaft  to  simulate  an  inviscid  flow,  should  give  the  result  that  the 
effective  velocity  everywhere  on  the  blade  is  1.0.  Previous  results  for  this  test  case 
[30]  ,  showed  that  this  result  was  not  attained. 

The  improved  algorithm  for  the  calculation  of  induced  velocities  presented  in  this 
thesis  shows  large  improvements  over  previous  results.  Figure  6-4  shows  contours  of 
effective  velocity  for  the  new  (analytic)  and  old  (discrete)  cases.  The  velocity  should 
be  1.0  everywhere  on  the  blade.  The  analytic  algorithm  presented  in  this  thesis 
shows  better  results  than  the  discrete  scheme  previously  used.  These  improvements, 
while  not  perfect,  carry  through  into  the  calculation  of  the  propeller  performance  as 
compared  to  experiment,  shown  in  Table  6.1. 

A  final  way  to  view  the  positive  trend  is  to  examine  the  circulation  distribution 
on  the  blade  predicted  by  using  the  new  and  old  induced  velocity  algorithms  with 
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PUF/UNCLE3D  PUF/UNCLE3D  Experiment 
discrete  analytic 

Kt  0.16412  0.15267  0.1455 

IOXq  0.30045  0.28191  0.2810 

Table  6.1:  Results  predicted  from  UNCLE  and  IFLOW  coupling  of  propeller  4119 
in  axial,  uniform  inflow.  The  discrete  column  is  the  previous  method  used  to  calcu¬ 
late  time-averaged  induced  velocities,  and  the  analytic  column  is  using  the  method 
presented  in  this  thesis. 


Figure  6-4:  These  two  plots  show  the  axial  component  of  the  effective  velocity  for 
Propeller  4119  in  uniform  inflow,  coupled  with  the  UNCLE  flow  solver. 
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Figure  6-5:  Blade  circulations  are  shown  for  the  old  and  new  coupling  algorithms, 
compared  to  the  stand-alone  calculation,  which  is  the  control  in  this  experiment. 

the  circulation  distribution  predicted  for  a  perfect  effective  inflow  of  uniform  axial 
velocity.  A  plot  of  the  predicted  blade  circulations  is  shown  in  Figure  6-5.  Notice 
how  much  closer  the  analytic  algorithm  recovers  the  true  solution.  Both  coupled 
algorithms  are  inaccurate  at  the  hub,  because  the  BANS  code  has  grown  a  shaft 
boundary  layer.  This  delta  at  the  hub  does  highlight  the  effect  of  not  properly 
accounting  for  a  thin  boundary  layer,  in  that  neglecting  the  presence  of  the  shaft 
boundary  layer  will  under  predict  the  true  circulation  at  the  hub. 

6.1.2  IFLOW-41  Coupling  with  Propeller  4119 

The  RANS  code  IFLOW-41  was  coupled  with  PUF  to  perform  several  analysis.  The 
grid  is  a  simple  0-type  grid,  with  48  cells  axially,  32  cells  radially,  and  32  cells 
azimuthally.  Because  the  purpose  of  this  study  wasn’t  to  capture  boundary  layer 
effects,  the  maximum  Y'^  is  allowed  to  be  19.  The  grid  is  shown  in  Figure  6-6  . 
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Figure  6-6:  Infinite  shaft  grid  used  in  IFL0W41  of  size  48x32x32. 


The  pressure  residual  history  for  two  different  blade  lattice  sizes  is  shown  in  Figure 
6-7.  This  figure  shows  that  the  coupling  with  IFLOW41  is  quickly  convergent.  A 
spike  in  the  pressure  residual  indicates  that  a  new  PUF  solution  for  body  forces  was 
injected  into  the  RANS  domain  at  that  iteration.  The  fact  that  the  convergence  is 
not  dependent  on  the  blade  lattice  size  indicates  that  the  body  force  coupling  routine, 
which  is  now  internal  to  PUF,  does  not  introduce  spurious  results. 


Propeller  Lattice  Convergence  Study 

A  propeller  lattice  convergence  was  conducted  to  find  the  minimum  necessary  pro¬ 
peller  lattice  size  required  to  produce  converged  answers.  The  results  are  presented  in 
Figure  6-8.  These  indicate  that  good  propeller  convergence  is  achieve  with  a  15X15 
blade  lattice  size  in  PUF.  Recall  that  the  real  key  here  is  that  the  number  of  discrete 
blade  positions  is  6.  This  is  an  order  of  magnitude  reduction  over  previous  techniques. 


IFL0W41  -  PDF  Coupled  History 


Figure  6-7:  Pressure  residual  history  for  IFLOW41  coupled  with  PUF.  The  solid  and 
dashed  line  are  for  two  different  blade  grid  densities  used  on  the  propeller  code. 
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Figure  6-8;  Increasing  the  propeller  lattice  size  when  coupled  with  IFLOW-41  pro¬ 
duces  results  that  are  closer  to  the  theoretical  results.  This  study  shows  that  use  of 
a  15X15  blade  lattice  is  sufficient  to  converge  a  simple  propeller  geometry,  such  as 
Propeller  4119. 
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6.2  Propeller  4119  at  10  Degree  Inclination 


A  preliminary  test  of  the  full  wake  alignment  and  higher-order  wake  panel  algorithms 
is  presented  here.  The  stability  of  the  new  method  is  shown  in  that  the  solutions 
converge,  and  produce  results  that  follow  intuition.  Since  there  are  no  experimental 
results  to  compare  with,  this  coupled  RANS  run  shows  that  the  new  propeller  program 
PUF-15  works  well  with  the  RANS  code  UNCLE-3D. 

Figure  6-9  highlights  new  physical  phenomena  which  are  modeled  when  a  realistic 
wake  model  is  introduced:  the  wake  detaches  from  the  lee  side  of  the  shaft,  while 
effects  of  radial  contraction  due  to  flow  acceleration  are  still  captured.  Conversely, 
on  the  windward  side  of  the  shaft,  the  modeling  is  not  so  accurate.  Realistically,  the 
wake  piles  up  on  itself,  while  the  model  here  allows  the  wake  to  pierce  the  shaft  if  the 
wake  stepping  algortihm  oversteps  the  innermost  streamline. 

A  pictorial  showing  the  propeller-induced  high  axial  velocities  being  convected 
downstream  at  an  angle  to  the  horizontal  is  shown  in  Figure  6-10. 

The  effects  of  the  wake  alignment  scheme  are  clearly  shown  in  the  next  set  of 
figures.  Figure  6-11  shows  the  magnitude  of  the  first  harmonic  of  Cp  at  the  0.5,  0.7, 
and  0.9  radii.  Figure  6-12  shows  the  magnitude  of  the  first  harmonic  of  Cp  at  the 
0.5,  0.7,  and  0.9  radii.  Figure  6-13  shows  the  phase  of  the  first  harmonic  of  Cp  at  the 
0.5,  0.7,  and  0.9  radii. 


6.3  Propeller  4679  at  7.5  Degree  Inclination 

DTMB  Propeller  4679  was  tested  at  the  DTMB  High  Speed  Basin  [10].  The  three 
bladed  propeller  was  driven  from  downstream  with  a  strut/pod  openwater  drive  sys¬ 
tem,  which  was  inclined  at  7.5  degrees  with  respect  to  the  horizontal  [11].  A  compu¬ 
tational  rendering  of  the  physical  model  is  shown  in  Figure  6-14. 
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Wake  Radial  Contraction 


Figure  6-9:  Using  the  flowfield  produced  by  Propeller  4119  on  a  shaft  inclined  at  10 
degrees,  and  coupled  with  the  UNCLE-3D  RANS  code. 


Contours  of  Swirl  Velocity 

Propeller  41 19  at  1 0  Degree  inclined  angle 


Figure  6-10:  Propeller  4119  at  10  Degree  inclined  angle  showing  the  convection  of 
swirl  downstream  of  the  propeller  swept  volume. 
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PUF:  First  Harmonic  Cp  Results  at  r/R=0.9 
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Figure  6-12:  Magnitude  of  the  First  Harmonic  of  Propeller  4119  at  10  Degree  inclined 


Figure  6-13:  Phase  of  the  First  Harmonic  of  Propeller  4119  at  10  Degree  inclined 
angle  Cp. 
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Figure  6-14:  Computational  model  of  the  DTMB  downstream  driven  propeller  shaft. 
The  propeller  sits  on  the  rounded  ball  near  the  nose. 

Propeller  DTMB  4679 

Inclined  Angle  7.5  degrees 

Advance  Coeff.  1.078 

Shaft  Configuration  Downstream  driven  podded  shaft 

6.3.1  Nominal  Flow  Convergence 

A  grid  convergence  study  was  conducted  on  three  grids  of  increasing  fineness  to  ensure 
the  solution  was  grid  independent.  Figure  6-15  shows  the  log  of  the  pressure  residual. 
Because  no  measurements  were  presented  for  the  pressure  coefficient  along  the  shaft 
or  boundary  layer  profiles,  it  is  not  possible  to  ascertain  the  validity  of  the  RANS 
model. 

It  is  interesting  to  examine  the  contours  of  axial  velocity  for  the  nominal  flow  case 
shown  in  Figure  6-16.  Note  the  high  velocities  over  the  ball  on  the  shaft  where  the 
propeller  sits.  The  large  radial  gradients  here  show  that  the  full  detail  of  the  shaft 
must  be  modeled  to  accurately  reproduce  the  empirical  propeller  measurements.. 
The  axial  velocity  contours  shown  learly  indicate  that  modeling  the  inflow  as  purely 
uniform,  axial  flow  will  produce  erroneous  computational  predictions. 
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0  100  200  300 

Multigrid  Cycle  Number 

Figure  6-15:  Log  of  pressure  residuals  for  nominal,  axial  flow  around  DTMB  down¬ 
stream  driven  shaft. 


Figure  6-16:  The  nominal  flow  around  the  DTMB  downstream  driven  test  shaft 


Figure  6-17:  Propeller  4679  blade  lattice  convergence  study.  The  20X20  lattice  shows 
nearly  the  same  results  as  the  15X15  lattice,  indicating  that  the  solution  has  achieved 
grid  independence. 

6.3.2  Propeller  Solution  Convergence 

A  propeller  convergence  test  was  conducted  by  varying  the  lattice  size  of  the  propeller 
while  holding  the  RANS  grid  constant.  The  resulting  propeller  circulation  distribution 
plots  are  shown  in  Figure  6-17.  Propeller  lattice  convergence  is  shown  at  a  20X20 
lattice.  The  expected  flattening  of  the  circulation  distribution  is  not  seen  near  the 
hub  due  to  the  axial  curvature  of  the  hub. 


6.3.3  Aligned  and  Non-Aligned  Wake  Results 

An  experiment  was  conducted  with  a  propeller  code  with  non-aligned  wakes  and  a 
propeller  code  with  aligned  wakes,  both  coupled  with  the  RANS  codes  for  Propeller 
4679  at  7.5  degree  inclined  angle.  This  test  quantifies  the  effects  upon  the  unsteady 
forces  of  aligning  the  wake.  The  coupled  results,  presented  in  Table  6.3.3  show  similar 
trends  presented  earlier  for  stand-alone  PUF.  The  large  difference  in  the  Z  component 
of  the  unsteady  force  is  expected  as  the  propeller  is  inclined  with  respect  to  the  Y 
axis,  thus  the  Z  directed  force  is  the  unsteady  component  of  the  side  force. 
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Component 

Wake  Aligned 

Non-aligned  Wake 

Fx 

0.01581 

0.01602 

Phase  Fx 

183 

159 

Fy 

0.00718 

0.00726 

Phase  Fy 

16 

36 

Fz 

0.01313 

0.01138 

Phase  Fz 

26 

8 

Table  6.2:  Predicted  1st  Harmonic  force  results  for  P4679  at  7.5  degree  inclined  angle 
for  aligned  and  non-aligned  wakes. 

6.3.4  Blade  Loading  Numerical  Comparison 

The  plots  of  the  blade  loading  for  each  discrete  propeller  position  are  shown  in  Figure 
6-18.  From  the  figure,  it  is  clear  that  there  are  magnitude  and  phase  differences 
between  the  two  propeller  programs. 

6.3.5  Experimental  Comparison 

The  results  from  the  experiments  conducted  at  DTMB  contain  Cp  data  at  constant 
radial  positions  along  the  propeller,  for  both  the  suction  and  pressure  side.  The  data 
are  both  amplitude  and  phase  information.  This  data  allows  a  key  comparison  to  be 
made  with  the  fully-aligned  wake  method,  since  results  presented  in  previous  sections 
have  shown  that  the  wake  alignments  procedure  alters  both  the  magnitude  and  phase 
of  the  resulting  blade  forces. 


Aligned  Wake 
Non-aligned  Wake 


Figure  6-18:  Propeller  4679  Steady  Cp  at  r/J?  =  0.5 


Propeller  4679  at  7.5  Degree  Inclined  Angle 
Steady  Cp  at  r/R=0.5 


Figure  6-19:  Propeller  4679  Steady  Cp  at  r/R  =  0.5 


96 


Propeller  4679  at  7.5  Degree  Inclined  Angle 
Steady  Cp  at  r/R=0.7 
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Figure  6-20:  Propeller  4679  Steady  Cp  at  r/R  =  0.7 


Propeller  4679  at  7.5  Degree  Inclined  Angle 
Steady  Cp  at  r/R=0.9 
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Figure  6-21:  Propeller  4679  Steady  Cp  at  r/R  =  0.9 


Propeller  4679  at  7.5  Degree  Inclined  Angle 
1st  Harmonic  Cp  at  r/R=0.9 


Figure  6-26:  Propeller  4679  first  harmonic  Cp  at  r/R  —  0.9 


First  Harmonic  Phase  Comparison 


Figure  6-27:  Propeller  4679  first  harmonic  phase  of  Cp  at  rfR  =  0.9 
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Chapter  7 


Conclusions 


This  thesis  has  focused  on  design  challenges  facing  propulsion  engineers.  By  im¬ 
proving  upon  the  state  of  the  art  techniques  in  vortex-lattice  lifting-surface  methods 
and  coupled  RANS  methods  as  specifically  applied  to  complex  propulsor  design  and 
analysis,  this  thesis  has  advanced  the  frontier  of  knowledge,  and  will  allow  practicing 
engineers  to  more  fully  explore  their  design  space. 

There  are  still  a  great  number  of  challenges  facing  the  industry  and  the  scientific 
research  community  which  supports  it.  The  largest  challenge  is  the  promise  of  more 
fully  populated  flow  simulation  codes  and  empirical  testing  techniques  to  enhance  our 
understanding  of  the  physics  of  the  problem,  and  then  apply  these  new  insights  into 
rigorous  design  methodologies.  This  thesis  has  shown  how  localized  applications  of 
physics-based  modeling  and  sleight  of  mind  paradigm  shifts  in  discretization  thoughts 
have  led  to  noticeable  gains  in  efficiency  and  accuracy.  And  as  far  as  pushing  the 
state  of  the  art,  this  is  a  practical  course  to  achieving  evolutionary  progress.  But 
as  basic  science  advances  upon  revolutionary  insights  based  on  ever  refined  length 
scales,  time  scales,  and  densities,  there  must  be  at  some  point  order  imposed  to  feed 
back  into  the  engineer  that  builds  a  product. 

There  are, also  ,  many  specific  challenges  facing  future  researchers  in  the  state  of 
the  art  as  envisioned  today, 

1.  Unsteady  drag  force  modeling:  While  the  state  of  the  art  is  pushing  for 
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steady  drag  modeling  on  a  vortex-lattice  lifting-surface  model,  the  unsteady 
blade  drag,  especially  in  the  highly  loaded  internal  propulsion  units  gaining 
popularity,  remains  an  open  question.  Experimental  evidence  shows  that  even 
the  steady  numerical  models  do  not  perform  as  well  as  was  hoped. 

2.  Thickness  modelingrThis  thesis  ignored  the  thickness  of  the  blade  (usually 
modeled  as  source  singularities)  upon  the  blade  solution.  This  was  purposeful 
in  that  if  the  unsteady  thickness  load  coupling  effects  are  modeled,  then  the 
unsteady  viscous  load  coupling  effects  must  too  be  modeled.  Ignoring  both  in 
tandem  is  hypothesized  to  create  no  undue  logical  errors. 

3.  Entropy  generation  in  RANS:  The  fact  that  blade  drag  appears  as  an 
entropy  source  in  the  RANS  code  could  be  modeled  through  another  tweaking 
of  the  right  hand  side  of  the  RANS  equation.  However,  the  effects  upon  firmly- 
established  turbulent  boundary  layer  models  would  certainly  lead  to  an  active 
debate. 


7.1  The  Future 

The  future  is  bright  for  today’s  engineers.  Ever  faster  computational  resources  with 
ever  expanding  memory  make  problems  solvable  today  that  were  untenable  twenty 
years  ago.  This  thesis  has  provided  a  few  new  tools  to  examine  some  of  the  more 
interesting  aspects  of  the  design  space  that  the  digital  and  hardware  revolution  brings 
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Appendix  A 


Validation  Techniques  for  PUF 
Type  Codes 


This  appendix  seeks  to  layout  validation  techniques  which  have  been  used  with  suc¬ 
cess  at  the  MIT  Marine  Hydrodynamics  Laboratory.  Along  with  validation  to  code 
improvments,  sections  are  also  provided  which  list  best  practice  trouble  shooting 
methods. 

Note  that  the  best  validation  is  with  experimental  data,  but  in  the  absence  of 
good  data  for  all  cases,  there  are  some  numerical  techniques  to  verify  the  accuracy  of 
a  coupled  technique. 

A.l  Coupled  Validation 

The  best  validation  between  a  vortex-lattice  lifting-surface  code  and  a  volume  solver 
(RANS,  Euler)  is  the  eflFective  inflow  equals  1.0  test. 

A. 1.1  Effective  Inflow  Equals  1.0  Test 

Running  the  volume  solver  in  inviscid  mode  turns  the  expensive  volume  solver  into  a 
Laplace  Solver. 

1.  Set  farfield  upstream  fluid  velocity  to  1.0 
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2.  Couple  PUF  and  the  volume  solver 


At  the  end  of  the  coupling,  the  effective  velocity  everywhere  on  the  blade  should  be 
1.0  since 


Veff  =  VranS  ~  y induced  (-A--!) 

and  there  is  no  vorticity  in  the  inflow  which  interacts  with  the  blade  vorticity  to  cause 
a  radial  redistribution  of  the  propeller  inflow  profile. 

This  procedure  also  results  in  a  completely  aligned  wake,  and  so  can  be  used  to 
validate  wake  self-alignment  techniques,  too. 

Plots  of  interest  would  show 

•  Contour  plot  of  the  effective  velocity  Vx  on  the  blade  surface 

•  Wake  lattice  pitch  at  various  spanwise  positions 

A. 1.2  Ve  Test 

Between  the  propeller  code  and  the  RANS  code,  the  tangential  induced  velocity  (or 
swirl  velocity)  must  be  the  same.  An  easy  check,  then,  is  to  extract  the  swirl  velocity 
from  the  RANS  domain  at  the  blade  trailing  edge,  and  compare  this  with  the  swirl 
velocity  output  by  the  propeller  code.  Differences  are  easily  attributable  to  force 
non-dimensionalization  and  scaling  issues,  and  should  be  immediately  addressed. 

A. 1.3  Blade  Isolation  Test 

A  good  test  to  test  the  coupling  suite  is  to  make  an  inflow  that  has  a  pie  slice  of  a 
lower  velocity.  Using  this  contrived  flowfield  as  an  input  to  the  propeller  code,  the 
output  induced  velocities  from  the  propeller  code  should  show  a  spike  in  the  same 
blade  position. 

Another  similary  worthwhile  test  is  to  couple  the  RANS  and  the  propeller  code 
where  the  RANS  inflow  has  harmonic  content.  Once  again,  an  overlay  of  the  induced 
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velocities  output  by  the  propeller  program  and  the  inflow,  should  show  contour  levels 
which  are  in  phase. 


A.  1.4  Coupled  Problems 

In  generating  a  new  coupling,  most  of  the  problems  lie  in  the  non-dimensionalization 
of  forces. 


A.  1.5  Coupling  Code  Non-Dimensionalizations 

Length  Scale 

The  length  scales  are  independently  scaled  between  PUF  and  the  volume  solver  (la¬ 
belled  as  RANS  in  the  following  equations).  The  ratio  of  the  length  scales  is  given  as 
A. 


LpUF 

Lrans 


(A.2) 


The  L  in  both  codes  is  taken  as  the  same  measure,  usually  of  propeller  diameter. 
Therefore,  the  volume  ratio  is  given  as 


Apt7F  _  A^Ap^jVg  _  ^3 
^RANS  ^RANS 


(A.3) 


Rotation  Speed 

By  standard  definition,  the  advance  ratio,  J,  is  given  as 


(A.4) 
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Within  PUF,  V  is  always  Vship  which  is  assumed  to  be  1.0  and  the  propeller  diameter 
D  is  2.0.  Thus,  the  rotation  speed  in  PUF  is  given  as 


‘IT’PUF  = 


1 


This  translates  into  an  given  as 


flpUF 


TT 

J 


(A.5) 


(A.6) 


Force 

Start  from  the  fact  that  by  similarity  considerations,  the  non-dimensional  thrust 
coefficient  must  match  between  the  two  codes. 


^TraNS  ^TpuF 


(A.7) 


By  definition,  then,  the  force  in  the  RANS  domain  can  be  shown  to  be 


Frans  =  KtPrans'<^%ansFrans 
=  KtPrans  ( 

=  j^KrPRANsVsD'RANS 

Substitute  in  the  definition  for  Kt  which  is  known  from  PUF. 


Frans 


Tpjjf  Prans^s^^anr 

PJJ  F 

—  ^Epuf,  PRANsYI^ 

^PUF 

=  ^TpuFPRANsVg 


(A.9) 
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If  in  the  RANS  code,  p  =  1.0  and  Vs  =  1.0  then 


Frans  _  J_ 
FpUF 


(A.IO) 


Force  Density 

Since  most  RANS  codes  use  a  force  potential  in  the  formulation  of  the  Navier-Stokes 
equations,  it  is  useful  to  see  how  a  force  density  in  the  PUF  coordinate  system  is 
transformed  to  a  force  density  in  the  RANS  coordinate  system. 


FRCDENSRANS  = 

Frans 


FRCDENSRANS 

FRCDENSPUF 


(A.11) 
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Appendix  B 


Potential  -  RANS  Coupling 


B.l  Effective  Inflow 


All  propeller  design  and  analysis  takes  place  in  the  effective  inflow.  The  calculation 
of  the  effective  inflow  is  the  starting  point  for  any  propeller  design.  An  accurate 
knowledge  of  the  effective  inflow  is  critical  to  the  success  of  the  propeller  design.  The 
effective  inflow  is  obscured  by  the  fact  that  the  vorticity  present  on  the  propeller  blades 
and  the  vorticity  within  the  boundary  layer  where  the  propeller  operates  interacts  in 
a  highly  non-linear  fashion. 

While  the  nominal  inflow  is  the  inflow  field  when  there  is  no  propeller  present,  it 
is  not  equal  to  the  effective  inflow.  This  is  due  to  the  effects  of  the  propeller  upon  the 
vorticity  in  the  inflow.  This  can  be  seen  from  the  definition  of  the  vorticity  vector. 

a;  =  VxV  (B.l) 


A  non-radially-uniform  meridional  inflow  suggests  the  presence  of  a  tangential  com¬ 
ponent  of  vorticity. 


^  _  dVr  dV:, 
dx  dr 


(B.2) 


The  tangential  vorticity  can  be  thought  of  as  a  ring  which  circumscribes  the  propeller 
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shaft  (Newman,  1977).  The  propeller  induced  velocity  causes  a  contraction  of  the 
tangential  vorticity  ring.  As  this  ring  contracts,  Kelvin’s  thereom  requires  that  the 
strength  of  the  vorticity  remain  constant.  Therefore,  from  equation  (B.2),  it  is  evident 
there  is  a  change  in  ^  to  counteract  the  change  in  If  there  were  no  tangential 
vorticity  present  in  the  inflow,  though,  the  nominal  and  effective  inflows  are  exactly 
equal.  This  fact  is  a  standard  technique  used  to  validate  the  any  coupling  presented 
between  a  potential- based  propeller  code  and  an  RANS  code. 


B.2  The  Coupled  Hull  Flow  Problem 

To  accurately  design  a  propeller  requires  an  accurate  knowledge  of  the  ship’s  resis¬ 
tance.  The  increased  fluid  velocity  along  the  stern  of  the  vessel  due  to  the  presence 
of  the  operating  propeller  causes  an  increase  in  drag  since  there  is  more  incomplete 
pressure  recovery.  Historically,  the  values  for  the  resistance  of  the  ship  with  and 
without  a  propeller  operating  are  related  through  the  thrust  deduction  coefficient. 

RT  =  il-t)T  (B.3) 

Modern  computation  techniques  which  model  the  interaction  of  the  propeller  and 
the  fluid  flow  around  the  hull  solve  both  the  effective  inflow  and  thrust  deduction 
problems  by  nearly  exactly  solving  for  the  fluid  flow  over  the  hull  with  the  propeller 
operating.  In  this  method,  the  entire  submarine  and  duct,  if  present,  is  modelled  in 
the  computer.  The  shear  flow  is  exactly  computed  along  the  body,  and  the  correct 
propeller  interaction  effects  are  also  modelled. 


B.3  Implementation  in  RANS  Formulation 

The  effects  of  the  propeller  are  introduced  into  the  RANS  domain  as  body  forces, 
which  appear  as  source  terms  on  the  right  hand  side  of  the  RANS  equation.  Equation 
B.4 


no 


The  body  force  term  Fi  must  be  properly  non-dimensionalized  for  the  RANS  code  at 
hand.  Usually,  the  non-dimensionalization  is, 


*  pU^/L 


(B.5) 


Where  fi  is  the  force  per  unit  volume  from  the  propeller  code. 
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Appendix  C 


Propeller  Numeric  Routines 

C.l  Romberg  Integration 

There  are  several  advantages  to  the  use  of  the  Romberg  adaptive  quadrature  scheme 
with  the  Richardson  extrapolation  [26]: 

•  Improved  accuracy.  The  truncation  error  associated  with  the  normal  trapezoidal 
scheme  is  0{h^)  while  with  the  Romberg  scheme  the  order  of  the  error  is  reduced 
to  0(^4). 

•  Calculation  reuse.  Previously  calculated  terms  are  used  in  the  extrapolation 
step  reducing  computational  expense. 

•  Controlled  accuracy.  The  stepping  scheme  which  steps  the  solution  to  finer 
levels  can  be  halted  once  a  desired  convergence  level  is  reached. 

•  Fast  convergence.  Because  the  Romberg  scheme  contains  an  extrapolated  error 
estimate,  fewer  steps  are  needed  to  achieve  convergence. 

•  The  Romberg  integration  scheme  has  been  validated  to  perform  very  well  when 
the  higher  derivatives  of  the  intergrand  are  large,  which  is  the  case  for  vortex 
influence  functions. 

The  basic  formula  of  the  trapezoidal  integration  is  given  in  Equation  C.l. 
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(C.l) 


7i,fc  =  -xifio)  +  f{b)  +  2y^^f{a  +  jh)] 

j=i 


where 


h  =2'=-'-! 

The  prediction  via  application  of  the  trapezoidal  scheme  via  successive  spacing  re¬ 
finement  can  be  accelarated  through  extrapolated  estimates  of  the  coarser  approx¬ 
imations.  Forming  these  successive  extrapolations  and  solution  into  a  matrix  gives 
the  elements  of  the  classic  Romberg  matrix, 


T,.,  =  (C.2) 

The  first  column  of  the  Romberg  matrix  contains  the  approximation  of  the  integral 
and  the  extrapolation  for  each  step  is  along  the  diagonal.  The  relative  truncation 
error  at  the  end  of  any  refinement  step  is  given  by, 


e 


'  Ti,k 


(C.3) 


C.2  Tri-linear  Interpolation 

The  art  of  data  interplation  is  a  careful  balance  of  the  computation  time,  accuracy, 
and  stability  of  the  method.  Many  schemese  were  investigated  before  the  tri-linear 
scheme  was  selected: 

•  tricubic  spline  interpolation  (too  unpredictable  in  high  gradient  regions) 

•  constant  strength  ( too  inaccurate  for  coarse  RANS  grids) 
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•  corner-weighted  block  interpolation  (inaccurate) 

In  the  PUF-RANS  coupled  scheme,  the  major  interpolations  involve  interpolating 
the  velocities  defined  on  cell  centers  to  the  PUF  control  points  (Neumann  boundary 
condition  points),  which  are  usually  not  coincidental.  Originally,  PUF  utilized  a 
three  dimensional  cubic  spline  interpolation.  This  scheme  proved  unstable  due  to  the 
presence  of  large,  higher  order,  gradients  in  the  flowfield,  especially  in  the  boundary 
layer  region.  This  thesis  implemented  a  tri-linear  interpolation  scheme. 

Given  the  position  and  velocities  on  the  eight  node  points  for  the  block  where  the 
field  point  is  located. 


Xi,  Vi,  Zi,  VXi,  Wi,  VZi  where  i  =  1..8 

The  linear  assumption  assumes  that  the  velocity  in  any  direction  at  any  field  point 
in  the  block  is  given  by, 

VX  {x,  y,  z)  =  Ax  +  By  +  Cz  +  Dxy  +  Exz  +  Fyz  +  Gxyz  H 
VY (x,  y,  z)  =  Ax  +  By  +  Cz  +  Dxy  +  Exz  V  Fyz  -\-  Gxyz  -\-  H 
VZ{x,  y,  z)  =  Ax  By  +  Cz  -\-  Dxy  -\-  Exz  +  Fyz  Gxyz  -I-  H 

where  the  weight  coefficients  A-H  differ  for  each  velocity  direction 

The  weight  coefficients  are  found  for  each  of  the  three  directions  independently  by 
solving  the  8x8  coeflScient  matrix.  Standard  LAPACK  LU  decpmposition  and  back 
substitution  routines  were  used  for  this  thesis. 


C.3  Elliptical  Streamline  Grids 

One  of  the  first  problems  that  is  encounter  by  the  propeller  program  is  to  find  a 
circumferential  mean  flowfield  from  the  RANS  total  flowfield,  so  that  the  propeller 
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blade  lattice  can  be  placed  along  these  streamlines.  Previous  schemes  took  the  cir¬ 
cumferential  mean  velocity  through  the  RANS  domain.  The  problem  came  that  in 
regions  of  high  velocity,  the  propeller  blade  growing  algorithm  would  take  too  large 
of  a  step  and  go  through  a  solid  boundary.  This  produced  propellers  that  extended 
into  hubs,  and  into  ducts. 

The  solution  which  this  thesis  incorporated  was  to  use  the  2-D  outer  boundary 
of  RANS  domain  as  the  edges  of  an  invisicid  flowfield,  and  solve  Laplace’s  equation 
internal  to  the  domain  to  give  beautiful  streamlines  which  follow  the  boundaries,  and 
will  keep  the  propeller  blades  in  the  fluid  domain. 

The  goal  of  any  grid  generation  routine  is  to  logically  organize  a  set  of  {x,  y,  z) 
points  within  the  fluid  domain  of  interest  so  that  the  conservation  laws  can  be  applied 
on  these  discrete  points. 

Within  the  realm  of  numerical  grid  generation  routines,  there  are  several  types 
of  grids.  The  generation  of  a  particular  grid  type  is  by  dictated  the  later  use  of  the 
grid.  A  particular  type  of  grid  is  only  suited  to  one,  or  a  few,  of  the  many  numerical 
implementations  of  the  governing  fluid  dynamic  equations.  On  a  broad  level,  there 
are  two  types  of  grids:  structured  and  unstructured. 

Another  classification  of  grid  generation  routines  is  between  those  that  are  fixed 
and  those  that  are  adaptive.  In  a  fixed  grid,  the  gridpoints  are  fixed.  In  an  adaptive 
scheme,  the  grid  points  are  altered  as  the  flow  solution  evolves  to  cluster  gridpoints 
in  regions  of  high  flow  gradients,  near  shocks  for  instance.  With  adaptive  grids,  the 
quality  of  the  flow  solution  is  increased  due  to  the  fineness  of  the  grid  in  the  regions 
that  require  it.  However,  this  increase  in  accuracy  comes  with  a  computational,  and 
programming,  cost. 

There  are  a  wide  variety  of  grid  generation  techniques  available  to  the  fluid  dy- 
namicist,  each  of  which  embodies  a  unique  set  of  trade-offs.  And  it  is  not  merely 
hearsay  that  the  quality  of  the  mesh  grid  greatly  affects  the  quality  of  the  numeri¬ 
cal  flow  solution.  This  paper  seeks  to  explain  the  rudimentary  differences,  and  offer 
simple  algorithms  for  their  implementations. 
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C.3.1  Structured  Grid  Generation 


The  goal  of  a  structured  grid  generation  routine  is  to  map  gridpoints  distributed  in 
the  X,  y  physical  fluid  domain  to  a  7/  computational  domain.  In  the  computation 
domain  =  At/  =  1  to  allow  a  flnite  dilference  discretization  of  the  differential  form 
of  the  flow  equations.  As  the  flnite  difference  equations  must  also  be  transformed  from 
the  physical  to  the  computational  domain,  the  grid  and  flow  solution  are  intimately 
coupled  by  the  metrics  of  the  transformation.  The  requirements  for  a  structured  grid 
are  that 

1.  There  be  a  one-to-one  mapping  such  that  gridlines  do  not  cross  each  other. 

Also,  for  numerical  accuracy  of  the  solution  to  the  flow  equations  we  desire  that, 

1.  The  distribution  of  gridpoints  is  smooth. 

2.  There  is  a  minimum  of  grid  line  skew. 

3.  Grid  line  orthogonality  or  near  orthogonality  is  maintained. 

4.  Grid  points  are  concentrated  in  regions  of  high  flow  gradients. 

C.3.2  Techniques  of  Structured  Grid  Generation 

There  are  three  historically-used  techniques  for  generating  structured  grids:  conformal 
mapping,  algebraic  mapping,  and  the  elliptic  partial  differential  equation. 

Conformal  Mapping  involves  the  use  of  complex  variables  to  map  the  physical  do¬ 
main  to  the  complex  domain.  While  the  metrics  are  guaranteed  to  be  good, 
conformal  mapping  is  applicable  only  in  two  dimensional  regions.  Even  in  two 
dimensional  regions,  very  few  geometries  are  known  for  which  conformal  map¬ 
ping  functions  can  be  generated.  Conformal  grids  are  no  longer  used  in  industry. 

Algebraic  Mapping  uses  an  algebraic  transformation  to  map  between  the  physical 
and  computational  domain.  These  routines  are  very  fast  and  easy  to  implement. 
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However,  control  of  grid  skewness  and  smoothness  is  difficult,  and  discontinuities 
at  the  boundary  are  deadly  when  later  solving  boundary  value  partial  differential 
equations. 

Elliptic  PDE  grids  are  generated  by  solving  Poisson’s  equation  with  a  specified 
boundary  geometry  and  either  Dirichlet  or  Neumann  boundary  conditions  (re¬ 
call  that  with  Dirichlet  boundary  conditions  the  points  on  the  boundary  are 
specified,  and  with  Neumann  boundary  conditions  the  gradient,  or  normal,  at 
the  boundary  is  specified).  PDE  grid  generation  routines  are  also  called  bound¬ 
ary  fitted,  elliptical,  or  Laplacian  grids.  These  routines  are  by  far  the  most 
flexible  and  widely  used  strucutred  grid  method. 

The  following  two  sections  will  discuss  the  theory  of  elliptic  grid  generation,  and 
how  algebraic  grids  found  a  proper  place. 

C.4  Poisson,  or  Elliptical,  Grid  Generation 

The  beauty  of  Laplace’s  equation  is  its  minimization  principle.  It  can  be  shown  that 
given  a  function  ^  which  satisfies  Laplace’s  equation, 

^xx  d"  ^yy  —  0  (C.4) 

This  function  ^  also  minimizes  the  quantity  7,  where  I  is  given  by 


7  =  f  (C.5) 

Jn 

Note  that  mathematically  expresses  the  length  of  the  grid  spacing  in  the  ^t) 
domain.  Hence,  a  solution  to  Laplace’s  equation  leads  to  a  minimization  of  grid 
density.  The  very  minimum  is  of  course  equal  spacing  ! 

So,  now  all  we  have  to  do  is  solve  Laplace’s  equation  in  the  computational  domain 
to  minimize  ^  and  77.  Laplace’s  equation  in  two  dimensions  for  the  (^,  rf)  coordinates 
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IS 


=0  (C.6) 

Vxx  "1“  Vyy  ^  (^•'^) 


In  equations  (C.6)  and  (C.7)  x  and  y  are  the  independent  variables  and  ^  and  y 
are  the  dependent  variables.  This  implies  that  we  are  given  a  distribution  of  gridpoints 
in  the  {x,  y)  domain  and  alter  the  corresponding  (^,  77)  gridpoints  until  equations  (C.6) 
and  (C.7)  are  satisfied.  In  reality,  though,  we  exactly  know  the  (^,  rj)  position  of  every 
grid  point  in  the  ^77  computational  domain.  So  the  trick  is  to  interchange  the  role  of 
the  independent  and  dependent  variables  in  equations  (C.6)  and  (C.7)  such  that  we 
are  now  solving  for  the  {x,  y)  position  of  every  (^,  77)  gridpoint. 

It  may  be  helpful  to  follow  the  algebraic  steps  necessary  to  derive  the  new  version 
of  Poisson’s  equation.  We  started  by  making  the  assertion  that 


^  =^(a:,y) 
77  =77(2;,  7/) 

and  that 


(C.8) 


Obviously  we  can  state  that  a  relationship  exists  for  the  inverse  relationship. 


(C.9) 


By  inspection  we  can  relate  the  two  transformation  matrices  and  assert  that 
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(C.IO) 


where  J  is  the  Jacobian  of  the  transformation  between  (^,97)  and  {x,y)  defined  as 


J  =  x^yr,  -  Xr,y^  (C.ll) 

Now  we  are  ready  to  perform  the  algebraic  steps  necessary  to  show  how  ^xx  equals 
a  function  where  x  and  y  are  differentiated  with  respect  to  ^  and  y.  We  start  with 
the  [1,1]  element  of  equation  (C.IO). 


use  the  chain  rule  of  differentiation 

■'•(j), *'•(*), 


When  this  process  is  carried  out  again  to  find  ^yy,rjxx,Vyy  ^-nd  after  several  alge¬ 
braic  groupings  are  made,  the  following  final  expressions  for  x  and  y  in  terms  of  ^ 
and  77  result. 


ax^^  -  2Px^jj  +  ^Xr,r,  =  0 


(C.12) 
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+  IVm  =  0 


(C.13) 


where: 

a  =  xl^yl 

13  =  X^Xr,  +  y^Vr, 

7  =  a:|  +  y| 

Equations  (C.12)  and  (C.13)  are  sometimes  referred  to  as  the  Thompson  equa¬ 
tions  after  Joe  Thompson  at  Mississippi  State  University  who  first  discovered  their 
properties. 

C.4.1  Solution  Techniques  for  the  Elliptical  Grid  Equations 

We  can  solve  equations  (C.12)  and  (C.13)  for  x  and  y,  respectively  using  standard 
iterative  numerical  methods  for  solving  elliptical  PDE’s  such  as  Gauss  Siedel,  SLOR, 
multigrid,  etc  . . .  The  non-linearities  can  be  assumed  away  be  treating  them  as  linear. 
A  simplified  solution  technique  follows  the  following  logic: 

1.  Provide  an  initial  guess  for  the  {x,y)j^k  values  of  the  gridpoints. 

2.  Calculate  {a,  7)j,fc  at  each  gridpoint. 

3.  Keeping  (a,  /?,  fixed,  perform  a  set  number  of  Gauss-Siedel  iterations  on 
the  X  coordinate. 

4.  Keeping  (a,  jd,  'y)j^k  fixed,  perform  a  set  number  of  Gauss-Siedel  iterations  on 
the  y  coordinate. 

5.  Update  boundary  conditions,  be  they  Dirichlet  or  Neumann. 

6.  Check  for  convergence  through  a  suitable  technique  -  for  instance,  calculate  the 
log  of  the  Z/2  norm  of  the  residual  solution  vector. 

7.  Repeat  steps  2-6  until  the  solution  is  converged  to  the  desired  level  of  conver¬ 
gence. 
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Note  that  the  non-linearities  present  in  equations  (C.12)  and  (C.13)  slow  the 
convergence  rate  for  elementary  iterative  techniques.  More  sophisticated  non-linear 
solution  methods  such  as  Newton-Rhapson  can  provide  a  faster  convergence. 

Although  no  forcing  functions  were  used  in  this  implementation,  it  would  be  an 
easy  manner  to  add  them  in  to  handle  more  exotic  domain  geometries. 

C.4.2  Forced  Functions 

Solving  equations  (C.12)  and  (C.13)  guarantees  a  smooth,  one  to  one  mapping.  How¬ 
ever,  an  orthogonal  grid  is  not  guaranteed  and  there  is  no  ability  to  cluster  gridpoints 
inthe  desired  region.  Also,  the  smoothness  of  the  Laplacian  solution  leads  to  pinched 
grids  near  convex  boundaries,  and  stretched  grids  near  concave  boundaries.  These 
deficiencies  can  be  overcome  through  the  use  of  forcing  functions  on  the  right  hand 
side  of  the  PDE’s,  leading  exactly  to  Poisson’s  equation.  Equations  (C.12)  and  (C.13) 
with  forcing  functions  are  given  by: 


"i*  ^Xj^  —  tJ  {^Px^  QXj^^ 


(C.14) 


ay((  -  Wvtrt  +  iVm  =  +  Qvtt) 


(C.15) 


The  forcing  functions  P  and  Q  allow  for  clustering  near  lines  of  constant  and 
r]i  and  near  individual  points  The  expressions  for  P  and  Q  are: 


Pi^,  V)  =  -Y1  (c.i6) 


Qitv)  =  -  E  aisign{ri  —  —^^bisign{r]  —  rii)e  Vi)'^  (C.17) 
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where 

Cj,  bi  are  amplification  factors 

ai,bi>  0  attracts  lines  to 
ai,bi  <0  repel  lines  from 

Ci,  di  are  decay  factors. 

The  amplification  factors  cause  gridlines  to  cluster  near  the  gridline  of  interest. 
The  decay  factors  set  the  degree  to  which  the  forcing  decays  as  you  get  away  from 
the  gridline  of  interest.  There  is  no  way  to  specify  a  priori  values  for  Uj,  bi,  Ci,  di. 

With  the  introduction  of  the  forcing  functions,  one  to  one  gridpoint  mapping  with 
non-overlapping  gridlines  between  the  ^r]  and  xy  planes  are  no  longer  guaranteed. 
The  grid  must  be  carefully  checked  before  it  is  used  in  a  flow  solver. 

Orthogonality  at  the  boundary  can  be  artificially  introduced  by  requiring  that  the 
grid  points  on  the  body,  which  is  itself  a  line  of  constant  rj  (i.e.,  rjo),  be  orthogonal 
to  the  gridpoints  above  the  body  surface  at  r]i.  The  Dirichlet  boundary  condition 
mathematically  can  be  expressed  as: 

x^Xn  +  y^yr,  =  0  (C.18) 

Obviously  this  condition  can  be  applied  at  any  other  boundary,  too.  Of  course, 
at  every  boundary,  only  one  type  of  boundary  condition  can  be  applied. 
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Appendix  D 


Propeller  Terminology 


D.l 

[B] 

D 

DTMB 

HGAP 

J 

Kt 

Kq 

NTSR 


n 

n 

P 

Q 

r 

R 


Abbreviations 

Propeller  Blade-on-Blade  influence  matrix 
Propeller  diameter  (  =  2.0  in  PUF  ) 

David  Taylor  Model  Basin  (now  Naval  Sea  Systems  Com¬ 
mand,  Carderock  Division) 

Hub  gap,  distance  between  the  hub  and  the  propeller, 

only  useful  for  tip  driven  propeller 

advance  coefficient,  =  ^ 

—  T 

—  _ 0. _ 

pn^D^ 

Number  of  Time  Steps  per  Revolution,  PUF  input  pa¬ 
rameter 

propeller  revolution  rate  in  rps 
Surface  normal  vector 
Pressure 
Propeller  torque 
Local  radius 

Propeller  max  radius  (  =  1.0  in  PUF  ) 
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Rt 

T 

T 

y® 


y© 

K// 

^induced 

^nominal 

Vtotal 


Vs 

Wt 

Wn 

z 

V 


p 

A 
— * 
r 

Tb 

Ts 

u 


Total  resistance 
Barehull  resistance 
Propeller  thrust 
Velocity  in  PUF 
Velocity  in  RANS 
Effective  Inflow  Velocity 
Induced  Inflow  Velocity 
Nominal  Inflow  Velocity 
Total  Inflow  Velocity 
Ship  velocity  (  =1.0  in  PUF  ) 

Propeller  wake  on  blade  influence  matrix 

Taylor  wake  fraction 

Nominal  volumetric  wake  fraction 

Number  of  propeller  blades 

Gradient  Operator 

density 

distance  scaling  factor  between  PUF  and  RANS 

Circulation 

Bound  circulation 

Shed  circulation,  dumped  into  the  wake 
Rotation  rate 
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D.2  Propeller  Terms 


Effective  Wake 


Nominal  Wake 

Total  Wake 


Non-physical  fluid  velocity  that  is  the  total  fluid  velocity 
in  the  region  of  the  propeller  minus  the  induced  veloci¬ 
ties  due  to  the  propeller.  Most  potential-based  propeller 
solvers  solve  the  propeller  problem  in  an  effective  wake. 
Fluid  flow  in  the  region  of  the  propeller  at  the  design 
vehicle  speed  without  the  presence  of  the  propeller 
Fluid  velocities  in  the  region  of  the  propeller  with  the 
effects  of  the  propeller. 


I 

I 
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